 |
PUMA
219
Portable University Model of the Atmosphere
|
|
PUMA
219
Portable University Model of the Atmosphere
|
00001 module pumamod 00002 00003 !*********************************************! 00004 ! Portable University Model of the Atmosphere ! 00005 !*********************************************! 00006 ! Version: 17.0 16-Feb-2011 ! 00007 !*********************************************! 00008 ! Klaus Fraedrich ! 00009 ! Frank Lunkeit - Edilbert Kirk ! 00010 ! Frank Sielmann - Torben Kunz ! 00011 ! Hartmut Borth ! 00012 !*********************************************! 00013 ! Meteorologisches Institut ! 00014 ! KlimaCampus - Universitaet Hamburg ! 00015 !*********************************************! 00016 ! http://www.mi.uni-hamburg.de/puma ! 00017 !*********************************************! 00018 00019 !**************************************************************! 00020 ! The number of processes for processing on parallel machines ! 00021 ! NLAT/2 must be dividable by <npro>. npro can be set by the ! 00022 ! option -n <npro> when calling the puma executable ! 00023 ! This option is only available if the code is compiled with ! 00024 ! an mpi compiler. ! 00025 !**************************************************************! 00026 integer :: npro = 1 00027 00028 !**************************************************************! 00029 ! The horizontal resolution of PUMA is set by defining the ! 00030 ! number of latitudes <nlev> with the 1st. command line ! 00031 ! parameter and the number of levels with the 2nd. command ! 00032 ! parameter. A typical call for T42 is: ! 00033 ! puma.x 64 10 ! 00034 ! which sets nlat=64 and nlev=10 ! 00035 !**************************************************************! 00036 integer :: nlat = 32 00037 00038 !example values: 32, 48, 64, 128, 192, 256, 512, 1024 00039 !truncation: T21, T31, T42, T85, T127, T170, T341, T682 00040 00041 integer :: nlev = 10 00042 00043 !*****************************************************! 00044 ! Grid related paramters, which are computed from the ! 00045 ! command line arguments <nlat> and <nlev> ! 00046 ! Preset values are for T21 (nlat=32) and nlev=10 ! 00047 ! ****************************************************! 00048 00049 integer :: nlem = 9 ! Levels - 1 00050 integer :: nlep = 11 ! Levels + 1 00051 integer :: nlsq = 100 ! Levels squared 00052 00053 integer :: nlon = 64 ! Longitudes = 2 * latitudes 00054 integer :: nlah = 16 ! Half of latitudes 00055 integer :: ntru = 21 ! (nlon-1) / 3 00056 integer :: ntp1 = 22 ! ntru + 1 00057 integer :: nzom = 44 ! Number of zonal modes 00058 integer :: nrsp = 506 ! (ntru+1) * (ntru+2) 00059 integer :: ncsp = 253 ! nrsp / 2 00060 integer :: nspp = 506 ! nodes per process 00061 integer :: nesp = 506 ! number of extended modes 00062 00063 integer :: nlpp = 32 ! Latitudes per process 00064 integer :: nhpp = 16 ! Half latitudes per process 00065 integer :: nhor = 2048 ! Horizontal part 00066 integer :: nugp = 2048 ! Horizontal total 00067 integer :: npgp = 1024 ! Horizontal total packed words 00068 00069 integer :: nud = 6 ! I/O unit for diagnostic output 00070 00071 !***********! 00072 ! filenames ! 00073 !***********! 00074 character (256) :: puma_namelist = "puma_namelist" 00075 character (256) :: puma_output = "puma_output" 00076 character (256) :: puma_diag = "puma_diag" 00077 character (256) :: puma_restart = "puma_restart" 00078 character (256) :: puma_status = "puma_status" 00079 character (256) :: efficiency_dat = "efficiency.dat" 00080 character (256) :: ppp_puma_txt = "ppp-puma.txt" 00081 character (256) :: puma_sp_init = "puma_sp_init" 00082 00083 ! ***************************************************************** 00084 ! * For multiruns the instance number is appended to the filename * 00085 ! * e.g.: puma_namelist_1 puma_diag_1 etc. for instance # 1 * 00086 ! ***************************************************************** 00087 00088 ! **************************************************************** 00089 ! * Don't touch the following parameter definitions ! * 00090 ! **************************************************************** 00091 integer, parameter :: PUMA = 0 ! Model ID 00092 integer, parameter :: PLASIM = 1 ! Model ID 00093 00094 parameter(NROOT = 0) ! Master node 00095 00096 parameter(PI = 3.141592653589793D0) ! Pi 00097 parameter(TWOPI = PI + PI) ! 2 Pi 00098 00099 parameter(AKAP_EARTH = 0.286 ) ! Kappa Earth 00100 parameter(AKAP_MARS = 0.2273) ! Kappa Mars 00101 parameter(ALR_EARTH = 0.0065) ! Lapse rate Earth 00102 parameter(ALR_MARS = 0.0025) ! Lapse rate Mars 00103 parameter(GA_EARTH = 9.81) ! Gravity Earth 00104 parameter(GA_MARS = 3.74) ! Gravity Mars 00105 parameter(GASCON_EARTH = 287.0) ! Gas constant for dry air on Earth 00106 parameter(GASCON_MARS = 188.9) ! Gas constant for dry air on Mars 00107 parameter(PSURF_EARTH = 101100.0) ! Mean Surface pressure [Pa] on Earth 00108 ! Trenberth 1981, J. Geoph. Res., Vol.86, 5238-5246 00109 parameter(PLARAD_EARTH = 6371000.0) ! Earth radius 00110 parameter(PLARAD_MARS = 3397000.0) ! Mars radius 00111 parameter(SID_DAY_EARTH= 86164.) ! Siderial day Earth 23h 56m 04s 00112 parameter(SID_DAY_MARS = 88642.) ! Siderial day Mars 24h 37m 22s 00113 00114 parameter(WW_EARTH = TWOPI/SID_DAY_EARTH) ! reciprocal of time scale 00115 ! on Earth [1/sec] 00116 parameter(WW_MARS = TWOPI/SID_DAY_MARS) ! reciprocal of time scale 00117 ! on Mars [1/sec] 00118 00119 parameter(CV_EARTH = PLARAD_EARTH * WW_EARTH) ! Velocity scale on Earth [m/s] 00120 parameter(CV_MARS = PLARAD_MARS * WW_MARS) ! Velocity scale on Mars [m/s] 00121 00122 parameter(CT_EARTH = CV_EARTH*CV_EARTH/GASCON_EARTH) !Temperature scale [K] 00123 ! on Earth 00124 parameter(CT_MARS = CV_MARS*CV_MARS/GASCON_MARS) !Temperature scale [K] 00125 ! on Mars 00126 00127 parameter(PNU = 0.02) ! Time filter 00128 parameter(PNU21 = 1.0 - 2.0*PNU) ! Time filter 2 00129 00130 ! ***************************************************************** 00131 ! * EZ: Factor to multiply the spherical harmonic Y_(1,0) to get * 00132 ! * the non-dimensional planetary vorticity 2 sin(phi). In PUMA * 00133 ! * Y_(1,0) = sqrt(3/2)*sin(phi) (normalization factor 1/sqrt(2)).* 00134 ! * The time scale must be given by Tscale = 1/Omega * 00135 ! ***************************************************************** 00136 parameter(EZ = 1.632993161855452D0) ! ez = 1 / sqrt(3/8) 00137 00138 00139 ! ************************************************************** 00140 ! * Planetary parameters & Scales * 00141 ! * ----------------------------- * 00142 ! * The Puma model is formulated in non-dimensional form with * 00143 ! * the planetary radius as length scale and the reciprocal of * 00144 ! * the planetary rotation rate as time scale. The temperature * 00145 ! * scale is given by the geopotential scale divided by the * 00146 ! * gas constant. * 00147 ! * For the time scale the length of the siderial day is used * 00148 ! * as basic unit * 00149 ! * The parameters are initialized for Earth settings. They * 00150 ! * may be modified by the namelist file <puma_namelist> * 00151 ! * * 00152 ! * The scales are derived internal quantities * 00153 ! ************************************************************** 00154 real :: sid_day = SID_DAY_EARTH ! Length of sideral day [sec] on Earth 00155 real :: plarad = PLARAD_EARTH ! Planetary radius [m] on Earth 00156 real :: gascon = GASCON_EARTH ! Dry air gas consant [J/K kg] on Earth 00157 real :: akap = AKAP_EARTH ! Kappa [] on Earth 00158 real :: alr = ALR_EARTH ! average lapse rate [K/km] on Earth 00159 real :: ga = GA_EARTH ! Gravity [m/sec*sec] on Earth 00160 real :: psurf = PSURF_EARTH ! Mean surface pressure for EARTH [Pa] 00161 00162 real :: ww = WW_EARTH ! reciprocal of time scale [1/sec] (Omega) 00163 real :: cv = CV_EARTH ! velocity scale [m/sec] on Earth 00164 real :: ct = CT_EARTH ! temperature scale [K] on Earth 00165 00166 ! ************************** 00167 ! * Global Integer Scalars * 00168 ! ************************** 00169 00170 logical :: lrestart = .false. ! Existing "puma_restart" sets to .true. 00171 logical :: lselect = .false. ! true: disable some zonal waves 00172 logical :: lspecsel = .false. ! true: disable some spectral modes 00173 00174 integer :: model = PUMA 00175 00176 integer :: kick = 1 ! kick > 0 initializes eddy generation 00177 integer :: nafter = 0 ! write data interval 0: controlled by nwpd 00178 integer :: nwpd = 1 ! number of writes per day 00179 integer :: ncoeff = 0 ! number of modes to print 00180 integer :: ndel = 6 ! ndel 00181 integer :: ndiag = 12 ! write diagnostics interval 00182 integer :: ngui = 0 ! activate Graphical User Interface 00183 integer :: nkits = 3 ! number of initial timesteps 00184 integer :: nlevt = 9 ! tropospheric levels (set_vertical_grid) 00185 integer :: noutput = 1 ! global switch for output on (1) or off (0) 00186 integer :: nwspini = 1 ! write sp_init after initialization 00187 integer :: nrun = 0 ! if (nstop == 0) nstop = nstep + nrun 00188 integer :: nstep1 = 0 ! start step (for cpu statistics) 00189 integer :: nstep = -1 ! current timestep step 0: 01-Jan-0001 00:00 00190 integer :: nstop = 0 ! finishing timestep 00191 integer :: ntspd = 0 ! number of timesteps per day 0 = auto 00192 integer :: mpstep = 0 ! minutes per step 0 = automatic 00193 integer :: ncu = 0 ! check unit (debug output) 00194 integer :: nwrioro = 1 ! controls output of orography 00195 integer :: nextout = 0 ! 1: extended output (entropy production) 00196 integer :: nruido = 0 ! 1: global constant, temporal noise 00197 ! 2: spatio-temporal noise 00198 ! 3: spatio-temporal equator symmetric 00199 integer :: nseedlen = 0 ! length of random seed (set by lib call) 00200 integer :: nmonths = 0 ! Simulation time (1 month = 30 days) 00201 integer :: nyears = 1 ! simulation time (1 year = 360 days) 00202 integer :: nsponge = 0 ! 1: Create sponge layer 00203 integer :: nhelsua = 0 ! 1: Set up Held & Suarez T_R field 00204 ! instead of original PUMA T_R field 00205 ! 2: Set up Held & Suarez T_R field 00206 ! instead of original PUMA T_R field 00207 ! AND use latitudinally varying 00208 ! heating timescale in PUMA (H&Z(94)), 00209 ! irrelevant for PumaPreProcessor (ppp) 00210 ! 3: Use latitudinally varying 00211 ! heating timescale in PUMA (H&Z(94)), 00212 ! irrelevant for PumaPreProcessor (ppp) 00213 integer :: ndiagp = 0 ! 0/1 switch for grid point diabatic heating 00214 integer :: nconv = 0 ! 0/1 switch for convecive heating 00215 integer :: nvg = 0 ! type of vertical grid 00216 ! 0 = linear 00217 ! 1 = Scinocca & Haynes 00218 ! 2 = Polvani & Kushner 00219 integer :: nenergy = 0 ! energy diagnostics (on/off 1/0) 00220 integer :: nentropy= 0 ! entropy diagnostics (on/off 1/0) 00221 integer :: ndheat = 0 ! energy recycling (on/off 1/0) 00222 00223 integer :: nradcv = 0 ! use two restoration fields 00224 00225 00226 00227 ! *********************** 00228 ! * Global Real Scalars * 00229 ! *********************** 00230 00231 real :: alpha = 1.0 ! Williams filter factor 00232 real :: alrs = 0.0 ! stratospheric lapse rate [K/m] 00233 real :: delt ! 2 pi / ntspd timestep interval 00234 real :: delt2 ! 2 * delt 00235 real :: dtep = 60.0 ! delta T equator <-> pole [K] 00236 real :: dtns = -70.0 ! delta T north <-> south [K] 00237 real :: dtrop = 12000.0 ! Tropopause height [m] 00238 real :: dttrp = 2.0 ! Tropopause smoothing [K] 00239 real :: dtzz = 10.0 ! delta(Theta)/H additional lapserate in 00240 ! Held & Suarez T_R field 00241 real :: orofac = 1.0 ! factor to scale the orograpy 00242 real :: plavor = EZ ! planetary vorticity 00243 real :: psmean = PSURF_EARTH ! Mean of Ps on Earth 00244 real :: rotspd = 1.0 ! rotation speed 1.0 = normal Earth rotation 00245 real :: sigmax = 6.0e-7 ! sigma for top half level 00246 real :: tdiss = 0.25 ! diffusion time scale [days] 00247 real :: tac = 360.0 ! length of annual cycle [days] (0 = no cycle) 00248 real :: pac = 0.0 ! phase of the annual cycle [days] 00249 real :: tgr = 288.0 ! Ground Temperature in mean profile [K] 00250 real :: dvdiff = 0.0 ! vertical diffusion coefficient [m2/s] 00251 ! ! dvdiff =0. means no vertical diffusion 00252 real :: disp = 0.0 ! noise dispersion 00253 real :: tauta = 40.0 ! heating timescale far from surface 00254 real :: tauts = 4.0 ! heating timescale close to surface 00255 real :: pspon = 50. ! apply sponge layer where p < pspon 00256 ! ! pressure [Pa] 00257 real :: sponk = 0.5 ! max. damping coefficient for sponge layer, 00258 ! ! unit: [1/day] 00259 00260 ! ************************** 00261 ! * Global Spectral Arrays * 00262 ! ************************** 00263 00264 real, allocatable :: sd(:,:) ! Spectral Divergence 00265 real, allocatable :: sdd(:,:) ! Difference between instances 00266 real, allocatable :: st(:,:) ! Spectral Temperature 00267 real, allocatable :: std(:,:) ! Difference between instances 00268 real, allocatable :: st1(:,:) ! Spectral Temperature at t-1 (for NEXTOUT == 1) 00269 real, allocatable :: st2(:,:) ! Spectral Temperature at t-2 (for NEXTOUT == 1) 00270 real, allocatable :: sz(:,:) ! Spectral Vorticity 00271 real, allocatable :: szd(:,:) ! Difference between instances 00272 real, allocatable :: sp(:) ! Spectral Pressure (ln Ps) 00273 real, allocatable :: spd(:) ! Difference between instances 00274 real, allocatable :: sq(:,:) ! For compatibility with PlaSim 00275 real, allocatable :: sp1(:) ! Spectral Pressure at t-1 (for NEXTOUT == 1) 00276 real, allocatable :: sp2(:) ! Spectral Pressure at t-2 (for NEXTOUT == 1) 00277 real, allocatable :: so(:) ! Spectral Orography 00278 real, allocatable :: sr1(:,:) ! Spectral Restoration Temperature 00279 real, allocatable :: sr2(:,:) ! Spectral Restoration Temperature 00280 00281 real, allocatable :: sdp(:,:) ! Spectral Divergence Partial 00282 real, allocatable :: stp(:,:) ! Spectral Temperature Partial 00283 real, allocatable :: szp(:,:) ! Spectral Vorticity Partial 00284 real, allocatable :: spp(:) ! Spectral Pressure Partial 00285 real, allocatable :: sop(:) ! Spectral Orography Partial 00286 real, allocatable :: srp1(:,:)! Spectral Restoration Partial 00287 real, allocatable :: srp2(:,:)! Spectral Restoration Partial 00288 00289 real, allocatable :: sdt(:,:) ! Spectral Divergence Tendency 00290 real, allocatable :: stt(:,:) ! Spectral Temperature Tendency 00291 real, allocatable :: szt(:,:) ! Spectral Vorticity Tendency 00292 real, allocatable :: spt(:) ! Spectral Pressure Tendency 00293 00294 real, allocatable :: sdm(:,:) ! Spectral Divergence Minus 00295 real, allocatable :: stm(:,:) ! Spectral Temperature Minus 00296 real, allocatable :: szm(:,:) ! Spectral Vorticity Minus 00297 real, allocatable :: spm(:) ! Spectral Pressure Minus 00298 00299 real, allocatable :: sak(:) ! Hyper diffusion 00300 real, allocatable :: srcn(:) ! 1.0 / (n * (n+1)) 00301 real, allocatable :: span(:) ! Pressure for diagnostics 00302 real, allocatable :: spnorm(:)! Factors for output normalization 00303 00304 integer, allocatable :: nindex(:) ! Holds wavenumber 00305 integer, allocatable :: nscatsp(:) ! Used for reduce_scatter op 00306 integer, allocatable :: nselzw(:) ! Enable/disable selected zonal waves 00307 integer, allocatable :: nselsp(:) ! Enable/disable slected spectral modes 00308 00309 ! *************************** 00310 ! * Global Gridpoint Arrays * 00311 ! *************************** 00312 00313 real, allocatable :: gd(:,:) ! Divergence 00314 real, allocatable :: gt(:,:) ! Temperature 00315 real, allocatable :: gz(:,:) ! Vorticity 00316 real, allocatable :: gu(:,:) ! u * cos(phi) 00317 real, allocatable :: gv(:,:) ! v * sin(phi) 00318 real, allocatable :: gp(:) ! Ln(Ps) 00319 real, allocatable :: gq(:,:) ! For compatibilty with PlaSim 00320 real, allocatable :: gfu(:,:) ! Term Fu in Primitive Equations 00321 real, allocatable :: gfv(:,:) ! Term Fv in Primitive Equations 00322 real, allocatable :: gut(:,:) ! Term u * T 00323 real, allocatable :: gvt(:,:) ! Term v * T 00324 real, allocatable :: gke(:,:) ! Kinetic energy u * u + v * v 00325 real, allocatable :: gpj(:) ! d(Ln(Ps)) / d(mu) 00326 real, allocatable :: rcsq(:) ! 1 / cos2(phi) 00327 real, allocatable :: ruido(:,:,:)! noise (nlon,nlat,nlev) 00328 real, allocatable :: ruidop(:,:) ! noise partial (nhor,nlev) 00329 real, allocatable :: gtdamp(:,:) ! 3D reciprocal damping times [1/sec] 00330 ! for relaxation in grid point space 00331 ! for radiative restoration temperature 00332 ! (e.g. for Held&Suarez) 00333 real, allocatable :: gr1(:,:) ! constant radiative restoration time scale 00334 real, allocatable :: gr2(:,:) ! variable radiative restoration time scale 00335 real, allocatable :: gtdampc(:,:)! the same as gtdamp, but for convective 00336 ! restoration temperature 00337 real, allocatable :: gr1c(:,:) ! constant convective restoration time scale 00338 real, allocatable :: gr2c(:,:) ! variable convective restoration time scale 00339 00340 ! ********************* 00341 ! * Diagnostic Arrays * 00342 ! ********************* 00343 00344 integer, allocatable :: ndil(:) ! Set diagnostics level 00345 00346 real, allocatable :: csu(:,:) ! Cross section u [m/s] 00347 real, allocatable :: csv(:,:) ! Cross section v [m/s] 00348 real, allocatable :: cst(:,:) ! Cross section T [Celsius] 00349 00350 real,allocatable :: denergy(:,:) ! energy diagnostics 00351 real,allocatable :: dentropy(:,:) ! entropy diagnostics 00352 00353 ! ******************* 00354 ! * Latitude Arrays * 00355 ! ******************* 00356 00357 character (3),allocatable :: chlat(:) ! label for latitudes 00358 real (kind=8),allocatable :: sid(:) ! sin(phi) 00359 real (kind=8),allocatable :: gwd(:) ! Gaussian weight (phi) 00360 real,allocatable :: csq(:) ! cos2(phi) 00361 real,allocatable :: rcs(:) ! 1/cos(phi) 00362 00363 ! **************** 00364 ! * Level Arrays * 00365 ! **************** 00366 00367 real, allocatable :: t0(:) ! reference temperature 00368 real, allocatable :: t0d(:) ! vertical t0 gradient 00369 real, allocatable :: taur(:) ! tau R [days] 00370 real, allocatable :: tauf(:) ! tau F [days] 00371 real, allocatable :: damp(:) ! 1.0 / (2 Pi * taur) 00372 real, allocatable :: fric(:) ! 1.0 / (2 Pi * tauf ) 00373 00374 real, allocatable :: bm1(:,:,:) 00375 real, allocatable :: dsigma(:) 00376 real, allocatable :: rdsig(:) 00377 real, allocatable :: sigma(:) ! full level sigma 00378 real, allocatable :: sigmh(:) ! half level sigma 00379 real, allocatable :: tkp(:) 00380 real, allocatable :: c(:,:) 00381 real, allocatable :: xlphi(:,:) ! matrix Lphi (g) 00382 real, allocatable :: xlt(:,:) ! matrix LT (tau) 00383 00384 ! ****************** 00385 ! * Parallel Stuff * 00386 ! ****************** 00387 00388 integer :: myworld = 0 ! MPI variable 00389 integer :: mpinfo = 0 ! MPI variable 00390 integer :: mypid = 0 ! My Process Id 00391 real :: tmstart = 0.0 ! CPU time at start 00392 real :: tmstop = 0.0 ! CPU time at stop 00393 character(80), allocatable :: ympname(:) ! Processor name 00394 00395 00396 ! ********************** 00397 ! * Multirun variables * 00398 ! ********************** 00399 00400 integer :: mrworld = 0 ! MPI communication 00401 integer :: mrinfo = 0 ! MPI info 00402 integer :: mrpid = -1 ! MPI instance id 00403 integer :: mrnum = 0 ! MPI number of instances 00404 integer :: mintru = 0 ! Lowest resolution of all instances 00405 integer :: mrdim = 0 ! Exchange dimension (min. NRSP) 00406 integer :: nsync = 0 ! Synchronization on or off 00407 integer, allocatable :: mrtru(:) ! Truncations of members 00408 00409 real :: syncstr = 0.0 ! Coupling strength (0 .. 1) 00410 real :: synctime = 0.0 ! Coupling time [days] 00411 00412 ! ****************************************** 00413 ! * GUI (Graphical User Interface for X11) * 00414 ! ****************************************** 00415 00416 parameter (NPARCS = 10) ! Number of GUI parameters 00417 integer :: nguidbg = 0 ! Flag for GUI debug output 00418 integer :: nshutdown = 0 ! Flag for shutdown request 00419 integer :: ndatim(6) = -1 ! Date & time array 00420 real(kind=4) :: parc(NPARCS) ! Values of GUI parameters 00421 real(kind=4) :: crap(NPARCS) ! Backup of parc(NPARCS) 00422 logical :: ldtep = .FALSE. ! DTEP changed by GUI 00423 logical :: ldtns = .FALSE. ! DTNS changed by GUI 00424 character(len=32) :: yplanet = "Earth" 00425 00426 ! *************** 00427 ! * Random seed * 00428 ! *************** 00429 00430 integer :: seed(8) = 0 ! settable in namelist 00431 integer, allocatable :: meed(:) ! machine dependent seed 00432 real :: ganext = 0.0! y part of gaussian noise 00433 00434 end module pumamod 00435 00436 !***************! 00437 ! MODULE RADMOD ! 00438 !***************! 00439 00440 module radmod ! Dummy declaration for compatibility 00441 use pumamod ! with PLASIM (needed in guimod) 00442 end module radmod 00443 00444 00445 ! ***************** ! 00446 ! * MODULE PPPMOD * ! 00447 ! ***************** ! 00448 00449 module prepmod 00450 integer :: num_ppp = 0 00451 integer :: nlat_ppp = 0 00452 integer :: nlev_ppp = 0 00453 00454 type ppp_type 00455 character (80) :: name ! name of variable or array 00456 logical :: isint ! .true. for integer 00457 integer :: n ! length of vector (1 for scalar) 00458 integer, pointer :: pint ! pointer to integer value or array 00459 real , pointer :: preal ! pointer to real value or array 00460 end type ppp_type 00461 00462 type(ppp_type) :: ppp_tab(30) 00463 00464 interface 00465 subroutine ppp_def_int(pname,nvar,ndim) 00466 character (*) :: pname 00467 integer, target :: nvar 00468 integer :: ndim 00469 end subroutine ppp_def_int 00470 subroutine ppp_def_real(pname,rvar,ndim) 00471 character (*) :: pname 00472 real , target :: rvar(*) 00473 integer :: ndim 00474 end subroutine ppp_def_real 00475 end interface 00476 00477 end module prepmod 00478 00479 00480 ! ********************* 00481 ! * PROGRAM PUMA_MAIN * 00482 ! ********************* 00483 00484 program puma_main 00485 use pumamod 00486 00487 ! *********** 00488 ! * History * 00489 ! *********** 00490 00491 ! 1972 - W. Bourke: 00492 ! An efficient one-level primitive equation spectral model 00493 ! Mon. Weath. Rev., 100, pp. 683-689 00494 00495 ! 1975 - B.J. Hoskins and A.J. Simmons: 00496 ! A multi-layer spectral model and the semi-implicit method 00497 ! Qart. J. R. Met. Soc., 101, pp. 637-655 00498 00499 ! 1993 - I.N. James and J.P. Dodd: 00500 ! A Simplified Global Circulation Model 00501 ! Users' Manual, Dept. of Meteorology, University of Reading 00502 00503 ! 1998 - Klaus Fraedrich, Edilbert Kirk, Frank Lunkeit 00504 ! Portable University Model of the Atmosphere 00505 ! DKRZ Technical Report No. 16 00506 00507 ! 2009 - PUMA Version 16.0 00508 ! http://www.mi.uni-hamburg.de/puma 00509 00510 ! ****************** 00511 ! * Recent Changes * 00512 ! ****************** 00513 00514 ! 10-Jun-2002 - Puma Workshop - Documentation of subroutine SPECTRAL 00515 ! 04-Jul-2002 - Frank Lunkeit - Annual cycle 00516 ! 08-Jul-2002 - Edilbert Kirk - Factor for rotation speed 00517 ! 25-Sep-2002 - Puma Workshop - Documentation of subroutine CALCGP 00518 ! 11-Nov-2002 - Edilbert Kirk - Add Orography to output file 00519 ! 26-Feb-2003 - Edilbert Kirk - Read preprocessed initial file 00520 ! 07-Sep-2004 - Edilbert Kirk - Graphical User Interface 00521 ! 23-Aug-2006 - Torben Kunz - Held & Suarez forcing 00522 ! 23-Aug-2006 - Torben Kunz - new spacing schemes of sigma levels 00523 ! 23-Aug-2006 - Edilbert Kirk - individual selection of zonal waves 00524 ! 23-Aug-2006 - Edilbert Kirk - optimized Legendre trasnformation module 00525 ! 19-Feb-2007 - Edilbert Kirk - new flexible restart I/O 00526 ! 15-Sep-2009 - Edilbert Kirk - static arrays replaced by allocatable 00527 ! 15-Sep-2009 - Frank Lunkeit - diagnostics for entropy production 00528 ! 27-Sep-2010 - Edilbert Kirk - cleaned up ruido routines 00529 00530 call mpstart 00531 call setfilenames 00532 call opendiag 00533 call read_resolution 00534 call resolution 00535 if (mrnum == 2) then 00536 call mrdimensions 00537 endif 00538 call allocate_arrays 00539 call prolog 00540 call master 00541 call epilog 00542 call guistop 00543 call mpstop 00544 stop 00545 end program puma_main 00546 00547 00548 ! *************************** 00549 ! * SUBROUTINE SETFILENAMES * 00550 ! *************************** 00551 00552 subroutine setfilenames 00553 use pumamod 00554 00555 character (3) :: mrext 00556 00557 if (mrpid < 0) return ! no multirun 00558 00559 write(mrext,'("_",i2.2)') mrpid 00560 00561 puma_namelist = trim(puma_namelist ) // mrext 00562 puma_output = trim(puma_output ) // mrext 00563 puma_diag = trim(puma_diag ) // mrext 00564 puma_restart = trim(puma_restart ) // mrext 00565 puma_status = trim(puma_status ) // mrext 00566 efficiency_dat = trim(efficiency_dat ) // mrext 00567 ppp_puma_txt = trim(ppp_puma_txt ) // mrext 00568 puma_sp_init = trim(puma_sp_init ) // mrext 00569 00570 return 00571 end 00572 00573 00574 ! *********************** 00575 ! * SUBROUTINE OPENDIAG * 00576 ! *********************** 00577 00578 subroutine opendiag 00579 use pumamod 00580 00581 if (mypid == NROOT) then 00582 open(nud,file=puma_diag) 00583 endif 00584 00585 return 00586 end 00587 00588 00589 ! ****************************** 00590 ! * SUBROUTINE ALLOCATE_ARRAYS * 00591 ! ****************************** 00592 00593 subroutine allocate_arrays 00594 use pumamod 00595 00596 allocate(sd(nesp,nlev)) ; sd(:,:) = 0.0 ! Spectral Divergence 00597 allocate(st(nesp,nlev)) ; st(:,:) = 0.0 ! Spectral Temperature 00598 allocate(sz(nesp,nlev)) ; sz(:,:) = 0.0 ! Spectral Vorticity 00599 allocate(sp(nesp)) ; sp(:) = 0.0 ! Spectral Pressure (ln Ps) 00600 allocate(so(nesp)) ; so(:) = 0.0 ! Spectral Orography 00601 allocate(sr1(nesp,nlev)) ; sr1(:,:) = 0.0 ! Spectral Restoration Temperature 00602 allocate(sr2(nesp,nlev)) ; sr2(:,:) = 0.0 ! Spectral Restoration Temperature 00603 allocate(sdp(nspp,nlev)) ; sdp(:,:) = 0.0 ! Spectral Divergence Partial 00604 allocate(stp(nspp,nlev)) ; stp(:,:) = 0.0 ! Spectral Temperature Partial 00605 allocate(szp(nspp,nlev)) ; szp(:,:) = 0.0 ! Spectral Vorticity Partial 00606 allocate(spp(nspp)) ; spp(:) = 0.0 ! Spectral Pressure Partial 00607 allocate(sop(nspp)) ; sop(:) = 0.0 ! Spectral Orography Partial 00608 allocate(srp1(nspp,nlev)) ; srp1(:,:)= 0.0 ! Spectral Restoration Partial 00609 allocate(srp2(nspp,nlev)) ; srp2(:,:)= 0.0 ! Spectral Restoration Partial 00610 allocate(sdt(nspp,nlev)) ; sdt(:,:) = 0.0 ! Spectral Divergence Tendency 00611 allocate(stt(nspp,nlev)) ; stt(:,:) = 0.0 ! Spectral Temperature Tendency 00612 allocate(szt(nspp,nlev)) ; szt(:,:) = 0.0 ! Spectral Vorticity Tendency 00613 allocate(spt(nspp)) ; spt(:) = 0.0 ! Spectral Pressure Tendency 00614 allocate(sdm(nspp,nlev)) ; sdm(:,:) = 0.0 ! Spectral Divergence Minus 00615 allocate(stm(nspp,nlev)) ; stm(:,:) = 0.0 ! Spectral Temperature Minus 00616 allocate(szm(nspp,nlev)) ; szm(:,:) = 0.0 ! Spectral Vorticity Minus 00617 allocate(spm(nspp)) ; spm(:) = 0.0 ! Spectral Pressure Minus 00618 allocate(sak(nesp)) ; sak(:) = 0.0 ! Hyper diffusion 00619 allocate(srcn(nesp)) ; srcn(:) = 0.0 ! 1.0 / (n * (n+1)) 00620 allocate(span(nesp)) ; span(:) = 0.0 ! Pressure for diagnostics 00621 allocate(spnorm(nesp)) ; spnorm(:)= 0.0 ! Factors for output normalization 00622 00623 allocate(nindex(nesp)) ; nindex(:) = ntru ! Holds wavenumber 00624 allocate(nscatsp(npro)) ; nscatsp(:) = nspp ! Used for reduce_scatter op 00625 allocate(nselzw(0:ntru)) ; nselzw(:) = 1 ! Enable selected zonal waves 00626 allocate(nselsp(ncsp)) ; nselsp(:) = 1 ! Enable slected spectral modes 00627 00628 allocate(gd(nhor,nlev)) ; gd(:,:) = 0.0 ! Divergence 00629 allocate(gt(nhor,nlev)) ; gt(:,:) = 0.0 ! Temperature 00630 allocate(gz(nhor,nlev)) ; gz(:,:) = 0.0 ! Vorticity 00631 allocate(gu(nhor,nlev)) ; gu(:,:) = 0.0 ! u * cos(phi) 00632 allocate(gv(nhor,nlev)) ; gv(:,:) = 0.0 ! v * sin(phi) 00633 allocate(gp(nhor)) ; gp(:) = 0.0 ! Ln(Ps) 00634 allocate(gfu(nhor,nlev)) ; gfu(:,:) = 0.0 ! Term Fu in Primitive Equations 00635 allocate(gfv(nhor,nlev)) ; gfv(:,:) = 0.0 ! Term Fv in Primitive Equations 00636 allocate(gut(nhor,nlev)) ; gut(:,:) = 0.0 ! Term u * T 00637 allocate(gvt(nhor,nlev)) ; gvt(:,:) = 0.0 ! Term v * T 00638 allocate(gke(nhor,nlev)) ; gke(:,:) = 0.0 ! Kinetic energy u * u + v * v 00639 allocate(gpj(nhor)) ; gpj(:) = 0.0 ! d(Ln(Ps)) / d(mu) 00640 00641 00642 allocate(rcsq(nhor)) ; rcsq(:) = 0.0 ! 1 / cos2(phi) 00643 00644 allocate(ndil(nlev)) ; ndil(:) = 0 00645 allocate(csu(nlat,nlev)) ; csu(:,:) = 0.0 00646 allocate(csv(nlat,nlev)) ; csv(:,:) = 0.0 00647 allocate(cst(nlat,nlev)) ; cst(:,:) = 0.0 00648 00649 allocate(chlat(nlat)) ; chlat(:) = ' ' 00650 allocate(sid(nlat)) ; sid(:) = 0.0 ! sin(phi) 00651 allocate(gwd(nlat)) ; gwd(:) = 0.0 ! Gaussian weight (phi) 00652 allocate(csq(nlat)) ; csq(:) = 0.0 ! cos2(phi) 00653 allocate(rcs(nlat)) ; rcs(:) = 0.0 ! 1/cos(phi) 00654 00655 allocate(t0(nlev)) ; t0(:) = 250.0 ! reference temperature 00656 allocate(t0d(nlev)) ; t0d(:) = 0.0 ! vertical t0 gradient 00657 allocate(taur(nlev)) ; taur(:) = 0.0 ! tau R [days] 00658 allocate(tauf(nlev)) ; tauf(:) = 0.0 ! tau F [days] 00659 allocate(damp(nlev)) ; damp(:) = 0.0 ! 1.0 / (2 Pi * taur) 00660 allocate(fric(nlev)) ; fric(:) = 0.0 ! 1.0 / (2 Pi * tauf ) 00661 allocate(dsigma(nlev)) ; dsigma(:) = 0.0 00662 allocate(rdsig(nlev)) ; rdsig(:) = 0.0 00663 allocate(sigma(nlev)) ; sigma(:) = 0.0 00664 allocate(sigmh(nlev)) ; sigmh(:) = 0.0 00665 allocate(tkp(nlev)) ; tkp(:) = 0.0 00666 allocate(c(nlev,nlev)) ; c(:,:) = 0.0 00667 allocate(xlphi(nlev,nlev)) ; xlphi(:,:) = 0.0 ! matrix Lphi (g) 00668 allocate(xlt(nlev,nlev)) ; xlt(:,:) = 0.0 ! matrix LT (tau) 00669 allocate(bm1(nlev,nlev,0:NTRU)) ; bm1(:,:,:) = 0.0 00670 00671 if (mrnum == 2) then 00672 allocate(sdd(nesp,nlev)) ; sdd(:,:) = 0.0 00673 allocate(std(nesp,nlev)) ; std(:,:) = 0.0 00674 allocate(szd(nesp,nlev)) ; szd(:,:) = 0.0 00675 allocate(spd(nesp )) ; spd(: ) = 0.0 00676 endif 00677 00678 return 00679 end subroutine allocate_arrays 00680 00681 00682 ! ================= 00683 ! SUBROUTINE PROLOG 00684 ! ================= 00685 00686 subroutine prolog 00687 use pumamod 00688 00689 character( 8) :: cpuma = 'PUMA-II ' 00690 character(80) :: pumaversion = '16.0 (27-Sep-2010)' 00691 real :: zsig(nlon*nlat) 00692 00693 if (mypid == NROOT) then 00694 call cpu_time(tmstart) 00695 write(nud,'(/," ****************************************************")') 00696 write(nud,'(" * PUMA ",a43," *")') trim(pumaversion) 00697 write(nud,'(" ****************************************************")') 00698 if (mrnum == 0) then 00699 write(nud,'(" * NTRU =",i4," NLEV =",i4," NLON = ",i4," NLAT =",i4," *")') & 00700 NTRU,NLEV,NLON,NLAT 00701 else 00702 do jpid = 1 , mrnum 00703 write(nud,'(" * PID =",i4," NTRU =",i4," NLEV = ",i4," *")') & 00704 jpid-1,mrtru(jpid),NLEV 00705 enddo 00706 endif 00707 write(nud,'(" ****************************************************")') 00708 if (NPRO > 1) then 00709 write(nud,'(/," ****************************************************")') 00710 do jpro = 1 , NPRO 00711 write(nud,'(" * CPU",i4,1x,a40," *")') jpro-1,ympname(jpro) 00712 enddo 00713 write(nud,'(" ****************************************************")') 00714 endif 00715 call restart_ini(lrestart,puma_restart) 00716 call inigau(NLAT,sid,gwd) 00717 call inilat 00718 call legpri 00719 call readnl 00720 call ppp_interface 00721 call initpm 00722 call initsi 00723 call altlat(csq,NLAT) ! csq -> alternating grid 00724 if (ngui > 0) call guistart 00725 if (nrun == 0 .and. nstop > 0) nrun = nstop-nstep 00726 if (nrun == 0) nrun = ntspd * (nyears * 360 + nmonths * 30) 00727 call initrandom ! set random seed 00728 endif ! (mypid == NROOT) 00729 00730 call mpbci(nruido) 00731 call initruido ! allocate ruido arrays 00732 00733 00734 if (nshutdown > 0) return ! If something went wrong in the init routines 00735 00736 ! *********************** 00737 ! * broadcast & scatter * 00738 ! *********************** 00739 00740 call mpscdn(sid,NHPP) ! real (kind=8) 00741 call mpscdn(gwd,NHPP) ! real (kind=8) 00742 call mpscrn(csq,NLPP) 00743 00744 do jlat = 1 , NLPP 00745 rcsq(1+(jlat-1)*NLON:jlat*NLON) = 1.0 / csq(jlat) 00746 enddo 00747 00748 ! broadcast integer 00749 00750 call mpbci(kick ) ! add noise for kick > 0 00751 call mpbci(nafter ) ! write data interval [steps] 00752 call mpbci(nwpd ) ! write data interval [writes per day] 00753 call mpbci(ncoeff ) ! number of modes to print 00754 call mpbci(ndel ) ! ndel 00755 call mpbci(noutput ) ! global output switch 00756 call mpbci(ndiag ) ! write diagnostics interval 00757 call mpbci(ngui ) ! GUI on (1) or off (0) 00758 call mpbci(nkits ) ! number of initial timesteps 00759 call mpbci(nlevt ) ! tropospheric levels 00760 call mpbci(nrun ) ! if (nstop == 0) nstop = nstep + nrun 00761 call mpbci(nstep ) ! current timestep 00762 call mpbci(nstop ) ! finishing timestep 00763 call mpbci(ntspd ) ! number of timesteps per day 00764 call mpbci(mpstep ) ! minutes per step 00765 call mpbci(nyears ) ! simulation time 00766 call mpbci(nmonths ) ! simulation time 00767 call mpbci(nextout ) ! write extended output 00768 call mpbci(nsponge) ! Switch for sponge layer 00769 call mpbci(nhelsua) ! Held & Suarez forcing 00770 call mpbci(ndiagp) ! 0/1 switch for new grid point diabatic heating 00771 call mpbci(nconv) ! 0/1 switch for convective heating 00772 call mpbci(nvg ) ! Type of vertical grid 00773 call mpbci(nenergy) ! energy diagnostics 00774 call mpbci(nentropy) ! entropy diagnostics 00775 call mpbci(ndheat) ! energy recycling 00776 call mpbci(nradcv) ! use two restoration fields 00777 00778 ! broadcast logical 00779 00780 call mpbcl(lrestart) ! true: read restart file, false: initial run 00781 call mpbcl(lselect ) ! true: disable some zonal waves 00782 call mpbcl(lspecsel) ! true: disable some spectral modes 00783 00784 ! broadcast real 00785 00786 call mpbcr(ww ) 00787 call mpbcr(v_scl ) 00788 call mpbcr(ct ) 00789 00790 call mpbcr(sid_day ) 00791 call mpbcr(plarad ) 00792 call mpbcr(gascon ) 00793 call mpbcr(akap ) 00794 call mpbcr(alr ) 00795 call mpbcr(ga ) 00796 call mpbcr(psurf ) 00797 call mpbcr(alpha ) ! Williams factor for time filter 00798 call mpbcr(dtep ) ! equator-pole temperature difference 00799 call mpbcr(dtns ) 00800 call mpbcr(dtrop ) 00801 call mpbcr(dttrp ) 00802 call mpbcr(tdiss ) 00803 call mpbcr(tac ) 00804 call mpbcr(pac ) 00805 call mpbcr(plavor ) 00806 call mpbcr(rotspd ) 00807 call mpbcr(sigmax ) ! sigma of top half level 00808 call mpbcr(tgr ) 00809 call mpbcr(dvdiff ) 00810 call mpbcr(disp ) 00811 call mpbcr(tauta ) 00812 call mpbcr(tauts ) 00813 call mpbcr(pspon ) 00814 call mpbcr(sponk ) 00815 00816 ! broadcast integer arrays 00817 00818 call mpbcin(ndil ,NLEV) 00819 call mpbcin(nselzw,NTP1) 00820 00821 ! broadcast real arrays 00822 00823 call mpbcrn(damp ,NLEV) 00824 call mpbcrn(dsigma,NLEV) 00825 call mpbcrn(fric ,NLEV) 00826 call mpbcrn(rdsig ,NLEV) 00827 call mpbcrn(taur ,NLEV) 00828 call mpbcrn(sigma ,NLEV) 00829 call mpbcrn(sigmh ,NLEV) 00830 call mpbcrn(t0 ,NLEV) 00831 call mpbcrn(t0d ,NLEV) 00832 call mpbcrn(tauf ,NLEV) 00833 call mpbcrn(tkp ,NLEV) 00834 00835 call mpbcrn(c ,NLSQ) 00836 call mpbcrn(xlphi ,NLSQ) 00837 call mpbcrn(xlt ,NLSQ) 00838 00839 ! scatter integer arrays 00840 00841 call mpscin(nindex,NSPP) 00842 call mpscrn(srcn ,NSPP) 00843 call mpscrn(sak ,NSPP) 00844 00845 call legini(nlat,nlpp,nesp,nlev,plavor,sid,gwd) 00846 00847 if (lrestart) then 00848 call read_atmos_restart 00849 if (mypid == NROOT) then 00850 if (kick > 10) call noise(kick-10) 00851 endif 00852 else 00853 call initfd 00854 endif 00855 00856 if (mypid == NROOT) then 00857 call printseed ! either namelist, clock initialized or from restart file 00858 endif 00859 00860 ! broadcast spectral arrays 00861 00862 call mpbcrn(sp,NESP) 00863 call mpbcrn(sd,NESP*NLEV) 00864 call mpbcrn(st,NESP*NLEV) 00865 call mpbcrn(sz,NESP*NLEV) 00866 00867 ! scatter spectral arrays 00868 00869 call mpscsp(sd,sdp,NLEV) 00870 call mpscsp(st,stp,NLEV) 00871 call mpscsp(sz,szp,NLEV) 00872 call mpscsp(sr1,srp1,NLEV) 00873 call mpscsp(sr2,srp2,NLEV) 00874 call mpscsp(sp,spp,1) 00875 call mpscsp(so,sop,1) 00876 00877 ! scatter gridpoint arrays 00878 00879 if (nruido > 0) call mpscgp(ruido,ruidop,NLEV) 00880 00881 ! 00882 ! initialize energy and entropy diagnostics 00883 ! 00884 if(nenergy > 0) then 00885 allocate(denergy(NHOR,9)) 00886 denergy(:,:)=0. 00887 endif 00888 if(nentropy > 0) then 00889 allocate(dentropy(NHOR,9)) 00890 dentropy(:,:)=0. 00891 endif 00892 if(ndheat > 1 .and. mypid == NROOT) then 00893 open(9,file=efficiency_dat,form='formatted') 00894 endif 00895 ! 00896 ! write first service record containing sigma coordinates 00897 ! 00898 if (mypid == NROOT) then 00899 if (noutput > 0) then 00900 istep = nstep 00901 if (istep > 0) istep = istep + nafter ! next write after restart 00902 open(40,file=puma_output,form='unformatted') 00903 call ntomin(istep,imin,ihour,iday,imonth,iyear) 00904 zsig(1:nlev) = sigmh(:) 00905 zsig(nlev+1:) = 0.0 00906 write(40) 333,0,iyear*10000+imonth*100+iday,0,nlon,nlat,nlev,ntru 00907 write(40) zsig 00908 endif ! (noutput > 0) 00909 endif ! (mypid == NROOT) 00910 return 00911 end subroutine prolog 00912 00913 00914 !===================! 00915 ! SUBROUTINE MASTER ! 00916 !================== ! 00917 00918 subroutine master 00919 use pumamod 00920 00921 if (nshutdown > 0) return ! if something went wrong in prolog already 00922 00923 ! *************************** 00924 ! * short initial timesteps * 00925 ! *************************** 00926 00927 ikits = nkits 00928 do jkits = 1 , ikits 00929 delt = (TWOPI/ntspd) / (2**nkits) 00930 delt2 = delt + delt 00931 call gridpoint 00932 call makebm 00933 call spectral 00934 nkits = nkits - 1 00935 enddo 00936 00937 delt = TWOPI/ntspd 00938 delt2 = delt + delt 00939 call makebm 00940 00941 nstep1 = nstep ! remember 1.st timestep 00942 00943 do jstep = 1 , nrun 00944 nstep = nstep + 1 00945 call ntomin(nstep,ndatim(5),ndatim(4),ndatim(3),ndatim(2),ndatim(1)) 00946 00947 ! ************************************************************ 00948 ! * calculation of non-linear quantities in grid point space * 00949 ! ************************************************************ 00950 00951 call gridpoint 00952 00953 if (mypid == NROOT) then 00954 if (mod(nstep,nafter) == 0 .and. noutput > 0) call outsp 00955 if (mod(nstep,ndiag ) == 0 .or. ngui > 0) call diag 00956 if (ncu > 0) call checkunit 00957 endif 00958 if (ngui > 0) call guistep_puma 00959 00960 ! ****************************** 00961 ! * adiabatic part of timestep * 00962 ! ****************************** 00963 00964 call spectral 00965 if (mod(nstep,nafter) == 0 .and. noutput > 0) call outgp 00966 if (nshutdown > 0) return 00967 enddo 00968 return 00969 end subroutine master 00970 00971 00972 ! ================= 00973 ! SUBROUTINE EPILOG 00974 ! ================= 00975 00976 subroutine epilog 00977 use pumamod 00978 real (kind=8) :: zut,zst 00979 integer (kind=8) :: imem,ipr,ipf,isw,idr,idw 00980 00981 if (mypid == NROOT) close(40) ! close output file 00982 00983 ! write restart file 00984 00985 if (mypid == NROOT) then 00986 call restart_prepare(puma_status) 00987 sp(1) = psmean ! save psmean 00988 call put_restart_integer('nstep' ,nstep ) 00989 call put_restart_integer('nlat' ,NLAT ) 00990 call put_restart_integer('nlon' ,NLON ) 00991 call put_restart_integer('nlev' ,NLEV ) 00992 call put_restart_integer('nrsp' ,NRSP ) 00993 00994 ! Save current random number generator seed 00995 00996 call random_seed(get=meed) 00997 call put_restart_array('seed',meed,nseedlen,1,1) 00998 call put_restart_array('ganext',ganext,1,1,1) 00999 01000 call put_restart_array('sz' ,sz ,NRSP,NESP,NLEV) 01001 call put_restart_array('sd' ,sd ,NRSP,NESP,NLEV) 01002 call put_restart_array('st' ,st ,NRSP,NESP,NLEV) 01003 call put_restart_array('sr1',sr1,NRSP,NESP,NLEV) 01004 call put_restart_array('sr2',sr2,NRSP,NESP,NLEV) 01005 call put_restart_array('sp' ,sp ,NRSP,NESP, 1) 01006 call put_restart_array('so' ,so ,NRSP,NESP, 1) 01007 if (nruido > 0) then 01008 call put_restart_array('ruido',ruido,nugp,nugp,nlev) 01009 endif 01010 endif 01011 01012 call mpputsp('szm',szm,NSPP,NLEV) 01013 call mpputsp('sdm',sdm,NSPP,NLEV) 01014 call mpputsp('stm',stm,NSPP,NLEV) 01015 call mpputsp('spm',spm,NSPP, 1) 01016 01017 ! write gridpoint arrays 01018 01019 if (allocated(gr1)) then 01020 call mpputgp('gr1',gr1,nhor,nlev) 01021 endif 01022 if (allocated(gr2)) then 01023 call mpputgp('gr2',gr2,nhor,nlev) 01024 endif 01025 if (allocated(gtdamp)) then 01026 call mpputgp('gtdamp',gtdamp,nhor,nlev) 01027 endif 01028 01029 if (allocated(gr1c)) then 01030 call mpputgp('gr1c',gr1c,nhor,nlev) 01031 endif 01032 if (allocated(gr2c)) then 01033 call mpputgp('gr2c',gr2c,nhor,nlev) 01034 endif 01035 if (allocated(gtdampc)) then 01036 call mpputgp('gtdampc',gtdampc,nhor,nlev) 01037 endif 01038 01039 if (mypid == NROOT) then 01040 ! Get resource stats from function resources in file pumax.c 01041 ires = nresources(zut,zst,imem,ipr,ipf,isw,idr,idw) 01042 call cpu_time(tmstop) 01043 tmrun = tmstop - tmstart 01044 if (nstep > nstep1) then 01045 zspy = tmrun * 360.0 * real(ntspd) / (nstep - nstep1) ! sec / siy 01046 zypd = (24.0 * 3600.0 / zspy) ! siy / day 01047 write(nud,'(/,"****************************************")') 01048 if (zut > 0.0) & 01049 write(nud, '("* User time : ", f10.3," sec *")') zut 01050 if (zst > 0.0) & 01051 write(nud, '("* System time : ", f10.3," sec *")') zst 01052 if (zut + zst > 0.0) tmrun = zut + zst 01053 write(nud, '("* Total CPU time : ", f10.3," sec *")') tmrun 01054 if (imem > 0) & 01055 write(nud, '("* Memory usage : ", f10.3," MB *")') imem * 0.000001 01056 if (ipr > 0) & 01057 write(nud, '("* Page reclaims : ", i6," pages *")') ipr 01058 if (ipf > 0) & 01059 write(nud, '("* Page faults : ", i6," pages *")') ipf 01060 if (isw > 0) & 01061 write(nud, '("* Page swaps : ", i6," pages *")') isw 01062 if (idr > 0) & 01063 write(nud, '("* Disk read : ", i6," blocks *")') idr 01064 if (idw > 0) & 01065 write(nud, '("* Disk write : ", i6," blocks *")') idw 01066 write(nud,'("****************************************")') 01067 if (zspy < 600.0) then 01068 write(nud,'("* Seconds per sim year: ",i6,9x,"*")') nint(zspy) 01069 else if (zspy < 900000.0) then 01070 write(nud,'("* Minutes per sim year ",i6,9x,"*")') nint(zspy/60.0) 01071 else 01072 write(nud,'("* Days per sim year: ",i6,5x,"*")') nint(zspy/86400.0) 01073 endif 01074 write(nud,'("* Sim years per day :",i7,9x,"*")') nint(zypd) 01075 write(nud,'("****************************************")') 01076 endif 01077 endif 01078 01079 return 01080 end subroutine epilog 01081 01082 ! ============================= 01083 ! SUBROUTINE READ_ATMOS_RESTART 01084 ! ============================= 01085 01086 subroutine read_atmos_restart 01087 use pumamod 01088 01089 integer :: k = 0 01090 01091 ! read scalars and full spectral arrays 01092 01093 if (mypid == NROOT) then 01094 call get_restart_integer('nstep',nstep) 01095 call get_restart_array('seed',meed,nseedlen,1,1) 01096 call get_restart_array('ganext',ganext,1,1,1) 01097 call get_restart_array('sz' ,sz ,NRSP,NESP,NLEV) 01098 call get_restart_array('sd' ,sd ,NRSP,NESP,NLEV) 01099 call get_restart_array('st' ,st ,NRSP,NESP,NLEV) 01100 call get_restart_array('sr1',sr1,NRSP,NESP,NLEV) 01101 call get_restart_array('sr2',sr2,NRSP,NESP,NLEV) 01102 call get_restart_array('sp' ,sp ,NRSP,NESP, 1) 01103 call get_restart_array('so' ,so ,NRSP,NESP, 1) 01104 if (nruido > 0) then 01105 call get_restart_array('ruido',ruido,nugp,nugp,nlev) 01106 endif 01107 psmean = sp(1) 01108 sp(1) = 0.0 01109 call random_seed(put=meed) 01110 endif 01111 01112 call mpbci(nstep) ! broadcast current timestep 01113 call mpbcr(psmean) ! broadcast mean surface pressure 01114 01115 ! read and scatter spectral arrays 01116 01117 call mpgetsp('szm',szm,NSPP,NLEV) 01118 call mpgetsp('sdm',sdm,NSPP,NLEV) 01119 call mpgetsp('stm',stm,NSPP,NLEV) 01120 call mpgetsp('spm',spm,NSPP, 1) 01121 01122 ! allocate, read and scatter gridpoint arrays 01123 01124 if (mypid == NROOT) call varseek('gr1',ktmp) 01125 call mpbci(ktmp) 01126 if (ktmp > 0) then 01127 allocate(gr1(nhor,nlev)) 01128 call mpgetgp('gr1',gr1,nhor,nlev) 01129 endif 01130 if (mypid == NROOT) call varseek('gr2',ktmp) 01131 call mpbci(ktmp) 01132 if (ktmp > 0) then 01133 allocate(gr2(nhor,nlev)) 01134 call mpgetgp('gr2',gr2,nhor,nlev) 01135 endif 01136 if (mypid == NROOT) call varseek('gtdamp',ktmp) 01137 call mpbci(ktmp) 01138 if (ktmp > 0) then 01139 allocate(gtdamp(nhor,nlev)) 01140 call mpgetgp('gtdamp',gtdamp,nhor,nlev) 01141 endif 01142 if (mypid == NROOT) call varseek('gr1c',ktmp) 01143 call mpbci(ktmp) 01144 if (ktmp > 0) then 01145 allocate(gr1c(nhor,nlev)) 01146 call mpgetgp('gr1c',gr1c,nhor,nlev) 01147 endif 01148 if (mypid == NROOT) call varseek('gr2c',ktmp) 01149 call mpbci(ktmp) 01150 if (ktmp > 0) then 01151 allocate(gr2c(nhor,nlev)) 01152 call mpgetgp('gr2c',gr2c,nhor,nlev) 01153 endif 01154 if (mypid == NROOT) call varseek('gtdampc',ktmp) 01155 call mpbci(ktmp) 01156 if (ktmp > 0) then 01157 allocate(gtdampc(nhor,nlev)) 01158 call mpgetgp('gtdampc',gtdampc,nhor,nlev) 01159 endif 01160 01161 return 01162 end subroutine read_atmos_restart 01163 01164 ! ================= 01165 ! SUBROUTINE INITFD 01166 ! ================= 01167 01168 subroutine initfd 01169 use pumamod 01170 01171 if (nkits < 1) nkits = 1 01172 01173 ! Look for start data and read them if there 01174 01175 call read_surf(129,so, 1,iread1) 01176 call read_surf(134,sp, 1,iread2) 01177 call read_surf(121,sr1,NLEV,iread3) 01178 call read_surf(122,sr2,NLEV,iread4) 01179 call read_vargp(123,NLEV,iread123) 01180 if (mypid == NROOT .and. iread123 == 0) then 01181 if (nhelsua > 1) then 01182 write(nud,*) "*** ERROR no *_surf_0123.sra file for Held&Suarez" 01183 stop 01184 endif 01185 endif 01186 01187 if (ndiagp > 0) then 01188 call read_vargp(121,NLEV,iread121) 01189 call read_vargp(122,NLEV,iread122) 01190 if (.not. allocated(gtdamp)) then 01191 call read_vargp(123,NLEV,iread123) 01192 endif 01193 if (mypid == NROOT) then 01194 if (iread121==0 .or. iread122==0 .or. iread123==0) then 01195 write(nud,*) "*** ERROR not all fields (121,122,123) for grid point heating found" 01196 stop 01197 endif 01198 endif 01199 endif 01200 01201 if (nconv > 0) then 01202 call read_vargp(124,NLEV,iread124) 01203 call read_vargp(125,NLEV,iread125) 01204 call read_vargp(126,NLEV,iread126) 01205 if (mypid == NROOT) then 01206 if (iread124==0 .or. iread125==0 .or. iread126==0) then 01207 write(nud,*) "*** ERROR not all fields (124,125,126) for convective heating found" 01208 stop 01209 endif 01210 endif 01211 endif 01212 01213 if (mypid == NROOT) then 01214 if (iread1==0 .or. iread2==0 .or. iread3==0 .or. iread4==0) then 01215 call setzt ! setup for aqua-planet 01216 else 01217 psmean = psurf * exp(spnorm(1) * sp(1)) 01218 sp(1) = 0.0 01219 so(:) = so(:) / (cv * cv) ! descale from [m2/s2] 01220 sr1(:,:) = sr1(:,:) / ct ! descale from [K] 01221 sr2(:,:) = sr2(:,:) / ct ! descale from [K] 01222 sr1(1,:) = sr1(1,:) - t0(:) * sqrt(2.0) ! subtract profile 01223 write(nud,'(a,f8.2,a)') ' Mean of Ps = ',0.01*psmean, '[hPa]' 01224 endif 01225 endif 01226 01227 ! Add initial noise if wanted 01228 01229 if (mypid == NROOT) then 01230 call printprofile 01231 if (kick > 10) then 01232 call noise(kick-10) 01233 else 01234 call noise(kick) 01235 endif 01236 endif ! (mypid == NROOT) 01237 01238 call mpscsp(sp,spm,1) 01239 if (mypid == NROOT) then 01240 st(1,:) = sr1(1,:) 01241 stm(1,:) = sr1(1,:) 01242 sz(3,:) = plavor 01243 szm(3,:) = plavor 01244 endif 01245 return 01246 end 01247 01248 01249 ! ========================== 01250 ! SUBROUTINE READ_RESOLUTION 01251 ! ========================== 01252 01253 subroutine read_resolution 01254 use pumamod 01255 01256 character (80) :: ylat 01257 character (80) :: ylev 01258 01259 if (mypid == NROOT) then 01260 call get_command_argument(1,ylat) 01261 call get_command_argument(2,ylev) 01262 read(ylat,*) nlat 01263 read(ylev,*) nlev 01264 endif 01265 01266 call mpbci(nlat) 01267 call mpbci(nlev) 01268 return 01269 end 01270 01271 01272 ! ===================== 01273 ! SUBROUTINE RESOLUTION 01274 ! ===================== 01275 01276 subroutine resolution 01277 use pumamod 01278 01279 nlem = nlev - 1 01280 nlep = nlev + 1 01281 nlsq = nlev * nlev 01282 01283 nlon = nlat + nlat ! Longitudes 01284 nlah = nlat / 2 01285 nlpp = nlat / npro 01286 nhpp = nlah / npro 01287 nhor = nlon * nlpp 01288 nugp = nlon * nlat 01289 npgp = nugp / 2 01290 01291 ntru = (nlon - 1) / 3 01292 ntp1 = ntru + 1 01293 nzom = ntp1 + ntp1 01294 nrsp = (ntru + 1) * (ntru + 2) 01295 ncsp = nrsp / 2 01296 nspp = (nrsp + npro - 1) / npro 01297 nesp = nspp * npro 01298 01299 return 01300 end 01301 01302 01303 ! ================= 01304 ! SUBROUTINE READNL 01305 ! ================= 01306 01307 subroutine readnl 01308 use pumamod 01309 01310 ! This workaround is necessaray, because allocatable arrays are 01311 ! not allowed in namelists for FORTRAN versions < F2003 01312 01313 integer, parameter :: MAXLEV = 100 01314 integer, parameter :: MAXSELZW = 42 01315 integer, parameter :: MAXSELSP = ((MAXSELZW+1) * (MAXSELZW+2)) / 2 01316 integer :: nselect(0:MAXSELZW) = 1 ! NSELECT can be used up tp T42 01317 integer :: nspecsel(MAXSELSP) = 1 ! Default setting: all modes active 01318 integer :: ndl(MAXLEV) = 0 ! Diagnostics off 01319 real :: restim(MAXLEV) = 0.0 ! Tau R 01320 real :: sigmah(MAXLEV) = 0.0 ! Half level sigma 01321 real :: t0k(MAXLEV) = 250.0 ! Reference temperature 01322 real :: tfrc(MAXLEV) = 0.0 ! Tau F 01323 01324 namelist /inp/ & 01325 akap , alpha , alr , alrs , disp , dtep & 01326 , dtns , dtrop , dttrp , dtzz , dvdiff & 01327 , ga , gascon & 01328 , kick , mpstep , nafter , ncoeff , nconv , ncu & 01329 , ndel , ndheat , ndiag , ndiagp , ndl , nenergy & 01330 , nentropy, nextout , ngui , nguidbg , nhelsua , nkits & 01331 , nlevt , nmonths , noutput , nradcv , nruido , nrun & 01332 , nselect , nspecsel, nsponge , nstep , nstop , nsync & 01333 , ntspd , nvg , nwpd , nwspini , nyears & 01334 , orofac , pac , plarad , pspon , psurf , restim & 01335 , rotspd , seed , sid_day , sigmah , sigmax , sponk & 01336 , syncstr , synctime, t0k & 01337 , tac , tauta , tauts , tdiss , tfrc , tgr 01338 01339 open(13,file=puma_namelist,iostat=ios) 01340 if (ios == 0) then 01341 read (13,inp) 01342 close(13) 01343 endif 01344 01345 !--- modify basic scales according to namelist 01346 ww = TWOPI/sid_day ! reciprocal of time scale 1/Omega 01347 cv = plarad*ww ! velocity scale (velocity at the equator) 01348 ct = cv*cv/gascon ! temperature scale from hydrostatic equation 01349 if (ntspd == 0) ntspd = (24 * nlat) / 32 ! automatic 01350 if (mpstep > 0) ntspd = 1440 / mpstep 01351 mpstep = 1440 / ntspd 01352 nafter = ntspd ! daily output 01353 if (nwpd > 0 .and. nwpd <= ntspd) then 01354 nafter = ntspd / nwpd 01355 endif 01356 if (ndiag < 1) ndiag = ntspd * 10 ! every 10th. day 01357 01358 if (synctime > 0.0) syncstr = 1.0 / (TWOPI * synctime) 01359 01360 write(nud,inp) 01361 01362 itru = ntru 01363 if (itru > MAXSELZW) itru = MAXSELZW 01364 icsp = ncsp 01365 if (icsp > MAXSELSP) icsp = MAXSELSP 01366 ilev = nlev 01367 if (ilev > MAXLEV) ilev = MAXLEV 01368 01369 nselzw(0:itru) = nselect(0:itru) ! Copy values to allocated array 01370 nselsp(1:icsp) = nspecsel(1:icsp) 01371 ndil(1:ilev) = ndl(1:ilev) 01372 taur(1:ilev) = restim(1:ilev) 01373 tauf(1:ilev) = tfrc(1:ilev) 01374 sigmh(1:ilev) = sigmah(1:ilev) 01375 t0(1:ilev) = t0k(1:ilev) 01376 01377 return 01378 end 01379 01380 01381 ! ====================== 01382 ! SUBROUTINE PPP_DEF_INT 01383 ! ====================== 01384 01385 subroutine ppp_def_int(pname,nvar,ndim) 01386 use prepmod 01387 01388 character (*) :: pname 01389 integer,target :: nvar 01390 01391 num_ppp = num_ppp + 1 01392 ppp_tab(num_ppp)%name = '[' // trim(pname) // ']' 01393 ppp_tab(num_ppp)%isint = .true. 01394 ppp_tab(num_ppp)%n = ndim 01395 ppp_tab(num_ppp)%pint => nvar 01396 ppp_tab(num_ppp)%preal => null() 01397 return 01398 end subroutine ppp_def_int 01399 01400 01401 ! ======================= 01402 ! SUBROUTINE PPP_DEF_REAL 01403 ! ======================= 01404 01405 subroutine ppp_def_real(pname,rvar,ndim) 01406 use prepmod 01407 character (*) :: pname 01408 real ,target :: rvar 01409 01410 num_ppp = num_ppp + 1 01411 ppp_tab(num_ppp)%name = '[' // trim(pname) // ']' 01412 ppp_tab(num_ppp)%isint = .false. 01413 ppp_tab(num_ppp)%n = ndim 01414 ppp_tab(num_ppp)%pint => null() 01415 ppp_tab(num_ppp)%preal => rvar 01416 return 01417 end subroutine ppp_def_real 01418 01419 01420 subroutine ppp_read_i(a,ndim,nread) 01421 integer :: a(ndim) 01422 integer :: n 01423 01424 nread = 0 01425 read (15,*) n 01426 if (n < 1 .or. n > ndim) return 01427 read (15,*) a(1:n) 01428 nread = n 01429 return 01430 end 01431 01432 01433 subroutine ppp_read_r(a,ndim,nread) 01434 real :: a(ndim) 01435 integer :: n 01436 01437 nread = 0 01438 read (15,*) n 01439 if (n < 1 .or. n > ndim) return 01440 read (15,*) a(1:n) 01441 nread = n 01442 return 01443 end 01444 01445 01446 ! ======================== 01447 ! SUBROUTINE PPP_INTERFACE 01448 ! ======================== 01449 01450 subroutine ppp_interface 01451 use pumamod 01452 use prepmod 01453 logical :: lexist 01454 integer :: iostat 01455 integer :: n 01456 integer :: ivar 01457 character (80) :: yname 01458 01459 inquire(file=ppp_puma_txt,exist=lexist) 01460 if (.not. lexist) return 01461 01462 call ppp_def_int('NLAT',nlat_ppp,1) 01463 call ppp_def_int('NLEV',nlev_ppp,1) 01464 01465 call ppp_def_real('SIGMH',sigmh,nlev) 01466 01467 write(nud,*) "*******************************" 01468 write(nud,*) "* Reading file <",trim(ppp_puma_txt),"> *" 01469 write(nud,*) "*******************************" 01470 open (15,file=ppp_puma_txt) 01471 read (15,'(A)',iostat=iostat) yname 01472 do while (trim(yname) /= '[END]' .and. iostat == 0) 01473 do j = 1 , num_ppp 01474 if (trim(yname) == ppp_tab(j)%name) then 01475 if (ppp_tab(j)%isint) then 01476 call ppp_read_i(ppp_tab(j)%pint,ppp_tab(j)%n,iread) 01477 if (iread == 0) then 01478 write(nud,*) "*** ERROR reading ",trim(yname)," from ",trim(ppp_puma_txt) 01479 stop 01480 else if (iread == 1) then 01481 write(nud,'("* ",A," = ",I10," *")') yname(1:15),ppp_tab(j)%pint 01482 else 01483 write(nud,'("* ",A," :",I5," items *")') yname(1:15),iread 01484 endif 01485 else 01486 call ppp_read_r(ppp_tab(j)%preal,ppp_tab(j)%n,iread) 01487 if (iread == 0) then 01488 write(nud,*) "*** ERROR reading ",trim(yname)," from ",trim(ppp_puma_txt) 01489 stop 01490 else if (iread == 1) then 01491 write(nud,'("* ",A," = ",G10.4," *")') yname(1:15),ppp_tab(j)%preal 01492 else 01493 write(nud,'("* ",A," :",I5," items *")') yname(1:15),iread 01494 endif 01495 endif 01496 exit 01497 endif 01498 enddo 01499 read (15,'(A)',iostat=iostat) yname 01500 enddo 01501 if (nlat_ppp /= 0 .and. nlat_ppp /= nlat) then 01502 write(nud,*) "*** ERROR *** ERROR *** ERROR *** ERROR ***" 01503 write(nud,*) "# of latitudes mismatch in preprocessor PPP and PUMA" 01504 write(nud,*) "NLAT in PPP : ",nlat_ppp," <",trim(ppp_puma_txt),">" 01505 write(nud,*) "NLAT in PUMA : ",nlat 01506 write(nud,*) "Aborting ..." 01507 stop 01508 endif 01509 if (nlev_ppp /= 0 .and. nlev_ppp /= nlev) then 01510 write(nud,*) "*** ERROR *** ERROR *** ERROR *** ERROR ***" 01511 write(nud,*) "# of levels mismatch in preprocessor PPP and PUMA" 01512 write(nud,*) "NLEV in PPP : ",nlev_ppp," <",trim(ppp_puma_txt),">" 01513 write(nud,*) "NLEV in PUMA : ",nlev 01514 write(nud,*) "Aborting ..." 01515 stop 01516 endif 01517 write(nud,*) "*******************************" 01518 01519 return 01520 end subroutine ppp_interface 01521 01522 01523 ! ============================= 01524 ! SUBROUTINE SELECT_ZONAL_WAVES 01525 ! ============================= 01526 01527 subroutine select_zonal_waves 01528 use pumamod 01529 01530 if (sum(nselzw(:)) /= NTP1) then ! some wavenumbers disabled 01531 lselect = .true. 01532 endif 01533 return 01534 end 01535 01536 ! ================================ 01537 ! SUBROUTINE SELECT_SPECTRAL_MODES 01538 ! ================================ 01539 01540 subroutine select_spectral_modes 01541 use pumamod 01542 01543 if (sum(nselsp(:)) /= NCSP) then ! some modes disabled 01544 lspecsel = .true. 01545 endif 01546 return 01547 end 01548 01549 ! ===================== 01550 ! * SET VERTICAL GRID * 01551 ! ===================== 01552 01553 subroutine set_vertical_grid 01554 01555 use pumamod 01556 01557 if (sigmh(NLEV) /= 0.0) return ! Already read in from namelist INP 01558 01559 if (nvg == 1) then ! Scinocca & Haynes sigma levels 01560 01561 if (nlevt >= NLEV) then ! Security check for 'nlevt' 01562 write(nud,*) '*** ERROR *** nlevt >= NLEV' 01563 write(nud,*) 'Number of levels (NLEV): ',NLEV 01564 write(nud,*) 'Number of tropospheric levels (nlevt): ',nlevt 01565 endif 01566 01567 ! troposphere: linear spacing in sigma 01568 ! stratosphere: linear spacing in log(sigma) 01569 ! after (see their Appendix): 01570 ! Scinocca, J. F. and P. H. Haynes (1998): Dynamical forcing of 01571 ! stratospheric planetary waves by tropospheric baroclinic eddies. 01572 ! J. Atmos. Sci., 55 (14), 2361-2392 01573 01574 ! Here, zsigtran is set to sigma at dtrop (tropopause height for 01575 ! construction of restoration temperature field). If tgr=288.15K, 01576 ! ALR=0.0065K/km and dtrop=11.km, then zsigtran=0.223 (=0.1 in 01577 ! Scinocca and Haynes (1998)). 01578 ! A smoothing of the transition between linear and logarithmic 01579 ! spacing, as noted in Scinocca and Haynes (1998), is not yet 01580 ! implemented. 01581 01582 zsigtran = (1. - alr * dtrop / tgr)**(ga/(gascon*alr)) 01583 zsigmin = 1. - (1. - zsigtran) / real(nlevt) 01584 01585 do jlev=1,NLEV 01586 if (jlev == 1) then 01587 sigmh(jlev) = SIGMAX 01588 elseif (jlev > 1 .and. jlev < NLEV - nlevt) then 01589 sigmh(jlev) = exp((log(SIGMAX) - log(zsigtran)) & 01590 & / real(NLEV - nlevt - 1) * real(NLEV - nlevt - jlev) 01591 + log(zsigtran)) 01592 elseif (jlev >= NLEV - nlevt .and. jlev < NLEV - 1) then 01593 sigmh(jlev) = (zsigtran - zsigmin) / real(nlevt - 1) 01594 * real(NLEV - 1 - jlev) + zsigmin 01595 elseif (jlev == NLEV - 1) then 01596 sigmh(jlev) = zsigmin 01597 elseif (jlev == NLEV) then 01598 sigmh(jlev) = 1. 01599 endif 01600 enddo 01601 return ! case nvg == 1 finished 01602 else if (nvg == 2) then ! Polvani & Kushner sigma levels 01603 inl = int(real(NLEV)/(1.0 - sigmax**(1.0/5.0))) 01604 do jlev=1,NLEV 01605 sigmh(jlev) = (real(jlev + inl - NLEV) / real(inl))**5 01606 enddo 01607 return 01608 01609 ! Default (nvg == 0) : equidistant sigma levels 01610 01611 else 01612 do jlev = 1 , NLEV 01613 sigmh(jlev) = real(jlev) / real(NLEV) 01614 enddo 01615 endif 01616 01617 return 01618 end 01619 01620 01621 ! ================= 01622 ! SUBROUTINE INITPM 01623 ! ================= 01624 01625 subroutine initpm 01626 use pumamod 01627 01628 real (kind=8) :: radea,zakk,zzakk 01629 real :: zsigb ! sigma_b for Held & Suarez frictional 01630 ! and heating timescales 01631 01632 radea = plarad ! Planet radius in high precision 01633 plavor = EZ * rotspd ! Planetary vorticity 01634 01635 ! ************************************************************* 01636 ! * carries out all initialisation of model prior to running. * 01637 ! * major sections identified with comments. * 01638 ! * this s/r sets the model parameters and all resolution * 01639 ! * dependent quantities. * 01640 ! ************************************************************* 01641 01642 if (lrestart) nkits=0 01643 01644 ! **************************************************** 01645 ! * Check for enabling / disabling zonal wavenumbers * 01646 ! **************************************************** 01647 01648 call select_zonal_waves 01649 if (npro == 1) call select_spectral_modes 01650 01651 ! ********************* 01652 ! * set vertical grid * 01653 ! ********************* 01654 01655 call set_vertical_grid 01656 01657 dsigma(1 ) = sigmh(1) 01658 dsigma(2:NLEV) = sigmh(2:NLEV) - sigmh(1:NLEM) 01659 01660 rdsig(:) = 0.5 / dsigma(:) 01661 01662 sigma(1 ) = 0.5 * sigmh(1) 01663 sigma(2:NLEV) = 0.5 * (sigmh(1:NLEM) + sigmh(2:NLEV)) 01664 01665 ! Initialize profile of tau R if not set in namelist 01666 01667 if (taur(NLEV) == 0.0) then 01668 do jlev = 1 , NLEV 01669 taur(jlev) = 158.0 / PI * atan(1.0 - sigma(jlev)) 01670 if (taur(jlev) > 30.0) taur(jlev) = 30.0 01671 enddo 01672 endif 01673 01674 ! Initialize profile of tau F if not set in namelist 01675 01676 if (tauf(NLEV) == 0.0) then 01677 do jlev = 1 , NLEV 01678 if (sigma(jlev) > 0.8) then 01679 tauf(jlev) = exp(10.0 * (1.0 - sigma(jlev))) / 2.718 01680 endif 01681 enddo 01682 endif 01683 01684 ! Compute 1.0 / (2 Pi * tau) for efficient use in calculations 01685 ! A day is 2 Pi in non dimensional units using omega as scaling 01686 01687 where (taur(:) > 0.0) 01688 damp(:) = 1.0 / (TWOPI * taur(:)) 01689 endwhere 01690 01691 where (tauf(:) > 0.0) 01692 fric(:) = 1.0 / (TWOPI * tauf(:)) 01693 endwhere 01694 01695 if (nsponge == 1) call sponge 01696 01697 01698 ! annual cycle period and phase in timesteps 01699 01700 if (tac > 0.0) tac = TWOPI / (ntspd * tac) 01701 pac = pac * ntspd 01702 01703 ! compute internal diffusion parameter 01704 01705 jdelh = ndel/2 01706 if (tdiss > 0.0) then 01707 zakk = ww*(radea**ndel)/(TWOPI*tdiss*((NTRU*(NTRU+1.))**jdelh)) 01708 else 01709 zakk = 0.0 01710 endif 01711 zzakk = zakk / (ww*(radea**ndel)) 01712 01713 ! set coefficients which depend on wavenumber 01714 01715 zrsq2 = 1.0 / sqrt(2.0) 01716 01717 jr =-1 01718 jw = 0 01719 do jm=0,NTRU 01720 do jn=jm,NTRU 01721 jr=jr+2 01722 ji=jr+1 01723 jw=jw+1 01724 nindex(jr)=jn 01725 nindex(ji)=jn 01726 spnorm(jr)=zrsq2 01727 spnorm(ji)=zrsq2 01728 zsq = jn * (jn+1) 01729 if (jn > 0) then 01730 srcn(jr) = 1.0 / zsq 01731 srcn(ji) = srcn(jr) 01732 endif 01733 sak(jr) = -zzakk * zsq**jdelh 01734 sak(ji) = sak(jr) 01735 enddo 01736 zrsq2=-zrsq2 01737 enddo 01738 01739 ! finally make temperatures dimensionless 01740 01741 dtns = dtns / ct 01742 dtep = dtep / ct 01743 ! dttrp = dttrp / ct 01744 t0(:) = t0(:) / ct 01745 01746 ! print out 01747 01748 write(nud,8120) 01749 write(nud,8000) 01750 write(nud,8010) NLEV 01751 write(nud,8020) NTRU 01752 write(nud,8030) NLAT 01753 write(nud,8040) NLON 01754 if (zakk == 0.0) then 01755 write(nud,8060) 01756 else 01757 write(nud,8070) ndel 01758 write(nud,8080) 01759 write(nud,8090) zakk,ndel 01760 write(nud,8100) tdiss 01761 endif 01762 write(nud,8110) PNU 01763 write(nud,8000) 01764 write(nud,8120) 01765 return 01766 01767 8000 format('*****************************************************') 01768 8010 format('* NLEV = ',i6,' Number of levels *') 01769 8020 format('* NTRU = ',i6,' Triangular truncation *') 01770 8030 format('* NLAT = ',i6,' Number of latitudes *') 01771 8040 format('* NLON = ',i6,' Number of longitues *') 01772 8060 format('* No lateral dissipation *') 01773 8070 format('* ndel = ',i6,' Lateral dissipation *') 01774 8080 format('* on vorticity, divergence and temperature *') 01775 8090 format('* with diffusion coefficient = ',e13.4,' m**',i1,'/s *') 01776 8100 format('* e-folding time for smallest scale is ',f7.3,' days *') 01777 8110 format('* Robert time filter with parameter PNU =',f8.3,' *') 01778 8120 format(/) 01779 end 01780 01781 01782 ! ================= 01783 ! SUBROUTINE MAKEBM 01784 ! ================= 01785 01786 subroutine makebm 01787 use pumamod 01788 01789 zdeltsq = delt * delt 01790 01791 do jlev1 = 1 , NLEV 01792 do jlev2 = 1 , NLEV 01793 zaq = zdeltsq * (t0(jlev1) * dsigma(jlev2)& 01794 & + dot_product(xlphi(:,jlev1),xlt(jlev2,:))) 01795 bm1(jlev2,jlev1,1:NTRU) = zaq 01796 enddo 01797 enddo 01798 01799 do jn=1,NTRU 01800 do jlev = 1 , NLEV 01801 bm1(jlev,jlev,jn) = bm1(jlev,jlev,jn) + 1.0 / (jn*(jn+1)) 01802 enddo 01803 call minvers(bm1(1,1,jn),NLEV) 01804 enddo 01805 return 01806 end 01807 01808 ! ================= 01809 ! SUBROUTINE INITSI 01810 ! ================= 01811 01812 subroutine initsi 01813 use pumamod 01814 01815 ! ********************************************** 01816 ! * Initialisation of the Semi Implicit scheme * 01817 ! ********************************************** 01818 01819 dimension zalp(NLEV),zh(NLEV) 01820 dimension ztautk(NLEV,NLEV) 01821 dimension ztaudt(NLEV,NLEV) 01822 01823 tkp(:) = akap * t0(:) 01824 t0d(1:NLEM) = t0(2:NLEV) - t0(1:NLEM) 01825 01826 zalp(2:NLEV) = log(sigmh(2:NLEV)) - log(sigmh(1:NLEM)) 01827 01828 xlphi(:,:) = 0.0 01829 xlphi(1,1) = 1.0 01830 do jlev = 2 , NLEV 01831 xlphi(jlev,jlev) = 1.0 - zalp(jlev)*sigmh(jlev-1)/dsigma(jlev) 01832 xlphi(jlev,1:jlev-1) = zalp(jlev) 01833 enddo 01834 01835 do jlev = 1 , NLEV 01836 c(jlev,:) = xlphi(:,jlev) * (dsigma(jlev) / dsigma(:)) 01837 enddo 01838 01839 ! *********************** tkp(i) = t0(i) * AKAP 01840 ! * matrix xlt - part 1 * 01841 ! *********************** 01842 01843 do jlev = 1 , NLEV 01844 ztautk(:,jlev) = tkp(jlev) * c(:,jlev) 01845 enddo 01846 01847 ! ********************* dsigma(i) = sigmh(i) - sigmh(i-1) 01848 ! * matrix xlt part 2 * rdsig (i) = 0.5 / dsigma(i) 01849 ! ********************* 01850 01851 ztaudt(1,1) = 0.5 * t0d(1) * (sigmh(1) - 1.0) 01852 ztaudt(2:NLEV,1) = 0.5 * t0d(1) * dsigma(2:NLEV) 01853 01854 do j= 2 , NLEV 01855 do i = 1 , j-1 01856 ztaudt(i,j) = dsigma(i) * rdsig(j) & 01857 * (t0d(j-1) * (sigmh(j-1)-1.0) + t0d(j) * (sigmh(j)-1.0)) 01858 enddo 01859 ztaudt(j,j) = 0.5 & 01860 * (t0d(j-1) * sigmh(j-1) + t0d(j) * (sigmh(j)-1.0)) 01861 do i = j+1 , NLEV 01862 ztaudt(i,j) = dsigma(i) * rdsig(j) & 01863 * (t0d(j-1) * sigmh(j-1) + t0d(j) * sigmh(j) ) 01864 enddo 01865 enddo 01866 01867 xlt(:,:) = ztautk(:,:) + ztaudt(:,:) 01868 01869 ! xlt finished 01870 01871 zfctr=0.001*cv*cv/ga 01872 do jlev=1,NLEV 01873 zh(jlev) = dot_product(xlphi(:,jlev),t0(:)) * zfctr 01874 enddo 01875 01876 ! ********************************** 01877 ! * write out vertical information * 01878 ! ********************************** 01879 01880 ilev = min(NLEV,5) 01881 write(nud,9001) 01882 write(nud,9002) 01883 write(nud,9003) 01884 write(nud,9002) 01885 do jlev=1,NLEV 01886 write(nud,9004) jlev,sigma(jlev),t0(jlev)*ct,zh(jlev) 01887 enddo 01888 write(nud,9002) 01889 write(nud,9001) 01890 01891 ! matrix c 01892 01893 write(nud,9012) 01894 write(nud,9013) 'c',(jlev,jlev=1,ilev) 01895 write(nud,9012) 01896 do jlev=1,NLEV 01897 write(nud,9014) jlev,(c(i,jlev),i=1,ilev) 01898 enddo 01899 write(nud,9012) 01900 write(nud,9001) 01901 01902 ! matrix xlphi 01903 01904 write(nud,9012) 01905 write(nud,9013) 'xlphi',(jlev,jlev=1,ilev) 01906 write(nud,9012) 01907 do jlev=1,NLEV 01908 write(nud,9014) jlev,(xlphi(i,jlev),i=1,ilev) 01909 enddo 01910 write(nud,9012) 01911 write(nud,9001) 01912 return 01913 9001 format(/) 01914 9002 format(33('*')) 01915 9003 format('* Lv * Sigma Basic-T Height *') 01916 9004 format('*',i3,' * ',3f8.3,' *') 01917 9012 format(69('*')) 01918 9013 format('* Lv * ',a5,i7,4i12,' *') 01919 9014 format('*',i3,' * ',5f12.8,' *') 01920 end 01921 01922 ! ===================== 01923 ! SUBROUTINE INITRANDOM 01924 ! ===================== 01925 01926 subroutine initrandom 01927 use pumamod 01928 integer :: i, clock 01929 01930 ! Set random number generator seed 01931 01932 call random_seed(size=nseedlen) 01933 allocate(meed(nseedlen)) 01934 01935 ! Take seed from namelist parameter 'SEED' ? 01936 01937 if (seed(1) /= 0) then 01938 meed(:) = 0 01939 i = nseedlen 01940 if (i > 8) i = 8 01941 meed(1:i) = seed(1:i) 01942 else 01943 call system_clock(count=clock) 01944 meed(:) = clock + 37 * (/(i,i=1,nseedlen)/) 01945 endif 01946 call random_seed(put=meed) 01947 return 01948 end 01949 01950 ! ==================== 01951 ! SUBROUTINE PRINTSEED 01952 ! ==================== 01953 01954 subroutine printseed 01955 use pumamod 01956 integer :: i 01957 01958 write (nud,9020) 01959 write (nud,9010) 01960 do i = 1 , nseedlen 01961 write (nud,9000) i,meed(i) 01962 enddo 01963 write (nud,9010) 01964 write (nud,9020) 01965 return 01966 9000 format('* seed(',i1,') = ',i10,' *') 01967 9010 format('************************') 01968 9020 format(/) 01969 end 01970 01971 ! ==================== 01972 ! SUBROUTINE INITRUIDO 01973 ! ==================== 01974 01975 subroutine initruido 01976 use pumamod 01977 if (nruido > 0) then 01978 allocate(ruido(nlon,nlat,nlev)) 01979 allocate(ruidop(nhor,nlev)) 01980 ruido = 77 01981 ruidop = 88 01982 endif 01983 return 01984 end 01985 01986 ! ==================== 01987 ! SUBROUTINE STEPRUIDO 01988 ! ==================== 01989 01990 subroutine stepruido 01991 use pumamod 01992 real :: zr 01993 integer :: need(8) 01994 01995 if (mypid == NROOT) then 01996 if (nruido == 1) then 01997 zr = disp*gasdev() 01998 ruido(:,:,:) = zr 01999 elseif (nruido == 2) then 02000 do jlev=1,NLEV 02001 do jlat=1,NLAT 02002 do jlon=1,NLON 02003 ruido(jlon,jlat,jlev) = disp*gasdev() 02004 enddo 02005 enddo 02006 enddo 02007 elseif (nruido == 3) then 02008 do jlev=1,NLEV 02009 do jlat=1,NLAT,2 02010 do jlon=1,NLON 02011 ruido(jlon,jlat ,jlev) = disp*gasdev() 02012 ruido(jlon,jlat+1,jlev) = ruido(jlon,jlat,jlev) 02013 enddo 02014 enddo 02015 enddo 02016 endif 02017 endif ! (mypid == NROOT) 02018 02019 call mpscgp(ruido,ruidop,NLEV) 02020 call random_seed(get=need) 02021 return 02022 end 02023 02024 ! ================== 02025 ! SUBROUTINE MINVERS 02026 ! ================== 02027 02028 subroutine minvers(a,n) 02029 dimension a(n,n),b(n,n),indx(n) 02030 02031 b = 0.0 02032 do j = 1 , n 02033 b(j,j) = 1.0 02034 enddo 02035 call ludcmp(a,n,indx) 02036 do j = 1 , n 02037 call lubksb(a,n,indx,b(1,j)) 02038 enddo 02039 a = b 02040 return 02041 end 02042 02043 ! ================= 02044 ! SUBROUTINE LUBKSB 02045 ! ================= 02046 02047 subroutine lubksb(a,n,indx,b) 02048 dimension a(n,n),b(n),indx(n) 02049 k = 0 02050 do i = 1 , n 02051 l = indx(i) 02052 sum = b(l) 02053 b(l) = b(i) 02054 if (k > 0) then 02055 do j = k , i-1 02056 sum = sum - a(i,j) * b(j) 02057 enddo 02058 else if (sum /= 0.0) then 02059 k = i 02060 endif 02061 b(i) = sum 02062 enddo 02063 02064 do i = n , 1 , -1 02065 sum = b(i) 02066 do j = i+1 , n 02067 sum = sum - a(i,j) * b(j) 02068 enddo 02069 b(i) = sum / a(i,i) 02070 enddo 02071 return 02072 end 02073 02074 ! ================= 02075 ! SUBROUTINE LUDCMP 02076 ! ================= 02077 02078 subroutine ludcmp(a,n,indx) 02079 dimension a(n,n),indx(n),vv(n) 02080 02081 d = 1.0 02082 vv = 1.0 / maxval(abs(a),2) 02083 02084 do 19 j = 1 , n 02085 do i = 2 , j-1 02086 a(i,j) = a(i,j) - dot_product(a(i,1:i-1),a(1:i-1,j)) 02087 enddo 02088 aamax = 0.0 02089 do i = j , n 02090 if (j > 1) & 02091 & a(i,j) = a(i,j) - dot_product(a(i,1:j-1),a(1:j-1,j)) 02092 dum = vv(i) * abs(a(i,j)) 02093 if (dum .ge. aamax) then 02094 imax = i 02095 aamax = dum 02096 endif 02097 enddo 02098 if (j .ne. imax) then 02099 do 17 k = 1 , n 02100 dum = a(imax,k) 02101 a(imax,k) = a(j,k) 02102 a(j,k) = dum 02103 17 continue 02104 d = -d 02105 vv(imax) = vv(j) 02106 endif 02107 indx(j) = imax 02108 if (a(j,j) == 0.0) a(j,j) = tiny(a(j,j)) 02109 if (j < n) a(j+1:n,j) = a(j+1:n,j) / a(j,j) 02110 19 continue 02111 return 02112 end 02113 02114 ! ============================= 02115 ! SUBROUTINE FILTER_ZONAL_WAVES 02116 ! ============================= 02117 02118 subroutine filter_zonal_waves(pfc) 02119 use pumamod 02120 dimension pfc(2,NLON/2,NLPP) 02121 02122 do jlat = 1 , NLPP 02123 pfc(1,1:NTP1,jlat) = pfc(1,1:NTP1,jlat) * nselzw(:) 02124 pfc(2,1:NTP1,jlat) = pfc(2,1:NTP1,jlat) * nselzw(:) 02125 enddo 02126 02127 return 02128 end 02129 02130 02131 ! ================================ 02132 ! SUBROUTINE FILTER_SPECTRAL_MODES 02133 ! ================================ 02134 02135 subroutine filter_spectral_modes 02136 use pumamod 02137 02138 j = 0 02139 k = -1 02140 do m = 0 , NTRU 02141 do n = m , NTRU 02142 k = k + 2 02143 j = j + 1 02144 if (nselsp(j) == 0) then 02145 spp(k:k+1 ) = 0.0 02146 sdp(k:k+1,:) = 0.0 02147 stp(k:k+1,:) = 0.0 02148 spt(k:k+1 ) = 0.0 02149 sdt(k:k+1,:) = 0.0 02150 stt(k:k+1,:) = 0.0 02151 spm(k:k+1 ) = 0.0 02152 sdm(k:k+1,:) = 0.0 02153 stm(k:k+1,:) = 0.0 02154 srp1(k:k+1,:) = 0.0 02155 srp2(k:k+1,:) = 0.0 02156 if (n < NTRU) then 02157 szp(k+2:k+3,:) = 0.0 02158 szt(k+2:k+3,:) = 0.0 02159 szm(k+2:k+3,:) = 0.0 02160 endif 02161 endif 02162 enddo 02163 enddo 02164 02165 return 02166 end 02167 02168 02169 ! ================ 02170 ! SUBROUTINE NOISE 02171 ! ================ 02172 02173 subroutine noise(kickval) 02174 use pumamod 02175 02176 ! kickval = -1 : read ln(ps) from puma_sp_init 02177 ! kickval = 0 : model runs zonally symmetric with no eddies 02178 ! kickval = 1 : add white noise to ln(Ps) asymmetric hemispheres 02179 ! kickval = 2 : add white noise to ln(Ps) symmetric to the equator 02180 ! kickval = 3 : force mode(1,2) of ln(Ps) allowing reproducable runs 02181 ! kickval = 4 : add white noise to symmetric zonal wavenumbers 7 of ln(Ps) 02182 02183 integer :: kickval 02184 integer :: jsp, jsp1, jn, jm 02185 integer :: jr, ji, ins 02186 real :: zr, zi, zscale, zrand 02187 02188 zscale = 0.000001 ! amplitude of noise 02189 zr = 0.01 ! kickval=3 value for mode(1,2) real 02190 zi = 0.005 ! kickval=3 value for mode(1,2) imag 02191 02192 select case (kickval) 02193 case (-1) 02194 open(71, file=puma_sp_init,form='unformatted',iostat=iostat) 02195 if (iostat /= 0) then 02196 write(nud,*) ' *** kick=-1: needs file <',trim(puma_sp_init),'> ***' 02197 stop 02198 endif 02199 read(71,iostat=iostat) sp(:) 02200 if (iostat /= 0) then 02201 write(nud,*) ' *** error reading file <',trim(puma_sp_init),'> ***' 02202 stop 02203 endif 02204 close(71) 02205 write(nud,*) 'initial ln(ps) field read from <',trim(puma_sp_init),'>' 02206 return 02207 case (0) ! do nothing 02208 case (1) 02209 jsp1=2*NTP1+1 02210 do jsp=jsp1,NRSP 02211 call random_number(zrand) 02212 if (mrpid > 0) zrand = zrand + mrpid * 0.01 02213 sp(jsp)=sp(jsp)+zscale*(zrand-0.5) 02214 enddo 02215 write(nud,*) 'white noise added' 02216 case (2) 02217 jr=2*NTP1-1 02218 do jm=1,NTRU 02219 do jn=jm,NTRU 02220 jr=jr+2 02221 ji=jr+1 02222 if (mod(jn+jm,2) == 0) then 02223 call random_number(zrand) 02224 if (mrpid > 0) zrand = zrand + mrpid * 0.01 02225 sp(jr)=sp(jr)+zscale*(zrand-0.5) 02226 sp(ji)=sp(ji)+zscale*(zrand-0.5) 02227 endif 02228 enddo 02229 enddo 02230 write(nud,*) 'symmetric white noise added' 02231 case (3) 02232 sp(2*NTP1+3) = sp(2*NTP1+3) + zr 02233 sp(2*NTP1+4) = sp(2*NTP1+4) + zi 02234 write(nud,*) 'mode(1,2) of ln(Ps) set to (',sp(2*NTP1+3),',',sp(2*NTP1+4),')' 02235 case (4) 02236 jr=2*NTP1-1 02237 do jm=1,NTRU 02238 do jn=jm,NTRU 02239 jr=jr+2 02240 ji=jr+1 02241 if (mod(jn+jm,2) == 0 .and. jm == 7) then 02242 call random_number(zrand) 02243 sp(jr)=sp(jr)+zscale*(zrand-0.5) 02244 sp(ji)=sp(ji)+zscale*(zrand-0.5) 02245 endif 02246 enddo 02247 enddo 02248 write(nud,*) 'symmetric zonal wavenumbers 7 of ln(Ps) perturbed', & 02249 & 'with white noise.' 02250 case default 02251 write(nud,*) 'Value ',kickval ,' for kickval not implemented.' 02252 stop 02253 end select 02254 02255 if (nwspini == 1) then 02256 open(71, file=puma_sp_init, form='unformatted') 02257 write(71) sp(:) 02258 close(71) 02259 endif 02260 02261 return 02262 end 02263 02264 ! ================ 02265 ! SUBROUTINE SETZT 02266 ! ================ 02267 subroutine setzt 02268 use pumamod 02269 02270 ! ************************************************************* 02271 ! * Set up the restoration temperature fields sr1 and sr2 * 02272 ! * for aqua planet conditions. * 02273 ! * The temperature at sigma = 1 is <tgr>, entered in kelvin. * 02274 ! * The lapse rate of ALR K/m is assumed under the tropopause * 02275 ! * and zero above. The tropopause is defined by <dtrop>. * 02276 ! * The smoothing ot the tropopause depends on <dttrp>. * 02277 ! ************************************************************* 02278 02279 dimension ztrs(NLEV) ! Mean profile 02280 dimension zfac(NLEV) 02281 02282 sr1(:,:) = 0.0 ! NESP,NLEV 02283 sr2(:,:) = 0.0 ! NESP,NLEV 02284 02285 ! Temperatures in [K] 02286 02287 zsigprev = 1.0 ! sigma value 02288 ztprev = tgr ! Temperature [K] 02289 zzprev = 0.0 ! Height [m] 02290 02291 do jlev = NLEV , 1 , -1 ! from bottom to top of atmosphere 02292 zzp=zzprev+(gascon*ztprev/ga)*log(zsigprev/sigma(jlev)) 02293 ztp=tgr-dtrop*alr ! temperature at tropopause 02294 ztp=ztp+sqrt((.5*alr*(zzp-dtrop))**2+dttrp**2) 02295 ztp=ztp-.5*alr*(zzp-dtrop) 02296 ztpm=.5*(ztprev+ztp) 02297 zzpp=zzprev+(gascon*ztpm/ga)*log(zsigprev/sigma(jlev)) 02298 ztpp=tgr-dtrop*alr 02299 ztpp=ztpp+sqrt((.5*alr*(zzpp-dtrop))**2+dttrp**2) 02300 ztpp=ztpp-.5*alr*(zzpp-dtrop) 02301 ztrs(jlev)=ztpp 02302 zzprev=zzprev+(.5*(ztpp+ztprev)*gascon/ga)*log(zsigprev/sigma(jlev)) 02303 ztprev=ztpp 02304 zsigprev=sigma(jlev) 02305 enddo 02306 02307 do jlev=1,NLEV 02308 ztrs(jlev)=ztrs(jlev)/ct 02309 enddo 02310 02311 !****************************************************************** 02312 ! loop to set array zfac - this controls temperature gradients as a 02313 ! function of sigma in tres. it is a sine wave from one at 02314 ! sigma = 1 to zero at stps (sigma at the tropopause) . 02315 !****************************************************************** 02316 ! first find sigma at dtrop 02317 ! 02318 zttrop=tgr-dtrop*alr 02319 ztps=(zttrop/tgr)**(ga/(alr*gascon)) 02320 ! 02321 ! now the latitudinal variation in tres is set up ( this being in terms 02322 ! of a deviation from t0 which is usually constant with height) 02323 ! 02324 zsqrt2 = sqrt(2.0) 02325 zsqrt04 = sqrt(0.4) 02326 zsqrt6 = sqrt(6.0) 02327 do 2100 jlev=1,NLEV 02328 zfac(jlev)=sin(0.5*PI*(sigma(jlev)-ztps)/(1.-ztps)) 02329 if (zfac(jlev).lt.0.0) zfac(jlev)=0.0 02330 sr1(1,jlev)=zsqrt2*(ztrs(jlev)-t0(jlev)) 02331 sr2(3,jlev)=(1./zsqrt6)*dtns*zfac(jlev) 02332 sr1(5,jlev)=-2./3.*zsqrt04*dtep*zfac(jlev) 02333 2100 continue 02334 write(nud,*) '**************************************************' 02335 write(nud,*) '* Restoration Temperature set up for aqua planet *' 02336 write(nud,*) '**************************************************' 02337 return 02338 end 02339 02340 ! ======================= 02341 ! SUBROUTINE PRINTPROFILE 02342 ! ======================= 02343 02344 subroutine printprofile 02345 use pumamod 02346 02347 ! ********************************** 02348 ! * write out vertical information * 02349 ! ********************************** 02350 02351 write(nud,9001) 02352 write(nud,9002) 02353 write(nud,9003) 02354 write(nud,9002) 02355 02356 do jlev=1,NLEV 02357 zt = (sr1(1,jlev)/sqrt(2.0) + t0(jlev)) * ct 02358 if (tauf(jlev) > 0.1) then 02359 write(nud,9004) jlev,sigma(jlev),zt,taur(jlev),tauf(jlev) 02360 else 02361 write(nud,9005) jlev,sigma(jlev),zt,taur(jlev) 02362 endif 02363 enddo 02364 02365 write(nud,9002) 02366 write(nud,9001) 02367 return 02368 9001 format(/) 02369 9002 format(36('*')) 02370 9003 format('* Lv * Sigma Restor-T tauR tauF *') 02371 9004 format('*',i3,' * ',f8.3,f9.3,2f5.1,' *') 02372 9005 format('*',i3,' * ',f8.3,f9.3,f5.1,' - *') 02373 end 02374 02375 02376 ! ==================== 02377 ! SUBROUTINE READ_SURF 02378 ! ==================== 02379 02380 subroutine read_surf(kcode,psp,klev,kread) 02381 use pumamod 02382 02383 logical :: lexist 02384 integer :: kread 02385 integer :: ihead(8) 02386 character(len=256) :: yfilename 02387 real :: psp(NESP,klev) 02388 real :: zgp(NUGP,klev) 02389 real :: zpp(NHOR,klev) 02390 02391 kread = 0 02392 if (mypid == NROOT) then 02393 if (NLAT < 1000) then 02394 write(yfilename,'("N",I3.3,"_surf_",I4.4,".sra")') NLAT,kcode 02395 else 02396 write(yfilename,'("N",I4.4,"_surf_",I4.4,".sra")') NLAT,kcode 02397 endif 02398 inquire(file=yfilename,exist=lexist) 02399 endif 02400 call mpbcl(lexist) 02401 if (.not. lexist) return 02402 02403 if (mypid == NROOT) then 02404 open(65,file=yfilename,form='formatted') 02405 write(nud,*) 'Reading file <',trim(yfilename),'>' 02406 do jlev = 1 , klev 02407 read (65,*) ihead(:) 02408 read (65,*) zgp(:,jlev) 02409 enddo 02410 close(65) 02411 if (kcode == 134) then 02412 write(nud,*) "Converting Ps to LnPs" 02413 zscale = log(100.0) - log(psurf) ! Input [hPa] / PSURF [Pa] 02414 zgp(:,:) = log(zgp(:,:)) + zscale 02415 endif 02416 call reg2alt(zgp,klev) 02417 endif ! (mypid == NROOT) 02418 call mpscgp(zgp,zpp,klev) 02419 call gp2fc(zpp,NLON,NLPP*klev) 02420 do jlev = 1 , klev 02421 call fc2sp(zpp(1,jlev),psp(1,jlev)) 02422 enddo 02423 call mpsum(psp,klev) 02424 kread = 1 02425 return 02426 end subroutine read_surf 02427 02428 02429 02430 ! ===================== 02431 ! SUBROUTINE READ_VARGP 02432 ! ===================== 02433 02434 subroutine read_vargp(kcode,klev,kread) 02435 use pumamod 02436 02437 logical :: lexist 02438 integer :: ihead(8) 02439 character(len=256) :: yfilename 02440 real :: zgp(NUGP,klev) 02441 02442 kread = 0 02443 if (mypid == NROOT) then 02444 if (NLAT < 1000) then 02445 write(yfilename,'("N",I3.3,"_surf_",I4.4,".sra")') NLAT,kcode 02446 else 02447 write(yfilename,'("N",I4.4,"_surf_",I4.4,".sra")') NLAT,kcode 02448 endif 02449 inquire(file=yfilename,exist=lexist) 02450 endif 02451 call mpbcl(lexist) 02452 if (.not. lexist) then 02453 if (mypid == NROOT) then 02454 write(nud,*) 'File <',trim(yfilename),'> not found' 02455 endif 02456 return 02457 endif 02458 02459 if (mypid == NROOT) then 02460 open(65,file=yfilename,form='formatted') 02461 write(nud,*) 'Reading file <',trim(yfilename),'>' 02462 do jlev = 1 , klev 02463 read (65,*) ihead(:) 02464 read (65,*) zgp(:,jlev) 02465 enddo 02466 close(65) 02467 call reg2alt(zgp,klev) 02468 endif ! (mypid == NROOT) 02469 02470 select case(kcode) 02471 case(121) 02472 !--- non-dimensionalize and shift const radiative rest. temp. 02473 if (mypid == NROOT) then 02474 zgp(:,:) = zgp(:,:)/ct 02475 do jhor = 1,nugp 02476 zgp(jhor,:) = zgp(jhor,:) - t0(:) 02477 enddo 02478 endif 02479 allocate(gr1(nhor,klev)) 02480 if (mypid == NROOT) then 02481 write(nud,*) 'Field gr1 allocated' 02482 endif 02483 call mpscgp(zgp,gr1,klev) 02484 case(122) 02485 !--- non-dimensionalize variable. radiative rest. temp. 02486 if (mypid == NROOT) then 02487 zgp(:,:) = zgp(:,:)/ct 02488 endif 02489 allocate(gr2(nhor,klev)) 02490 if (mypid == NROOT) then 02491 write(nud,*) 'Field gr2 allocated' 02492 endif 02493 call mpscgp(zgp,gr2,klev) 02494 case(123) 02495 !--- non-dimensionalize radiative relaxation time scale 02496 if (mypid == NROOT) then 02497 zgp(:,:) = zgp(:,:)/ww 02498 endif 02499 allocate(gtdamp(nhor,klev)) 02500 if (mypid == NROOT) then 02501 write(nud,*) 'Field gtdamp allocated' 02502 endif 02503 call mpscgp(zgp,gtdamp,klev) 02504 case(124) 02505 !--- non-dimensionalize and shift const. convective rest. temp. 02506 if (mypid == NROOT) then 02507 zgp(:,:) = zgp(:,:)/ct 02508 do jhor = 1,nugp 02509 zgp(jhor,:) = zgp(jhor,:) - t0(:) 02510 enddo 02511 endif 02512 allocate(gr1c(nhor,klev)) 02513 if (mypid == NROOT) then 02514 write(nud,*) 'Field gr1c allocated' 02515 endif 02516 call mpscgp(zgp,gr1c,klev) 02517 case(125) 02518 !--- non-dimensionalize variable. convective rest. temp. 02519 if (mypid == NROOT) then 02520 zgp(:,:) = zgp(:,:)/ct 02521 endif 02522 allocate(gr2c(nhor,klev)) 02523 if (mypid == NROOT) then 02524 write(nud,*) 'Field gr2c allocated' 02525 endif 02526 call mpscgp(zgp,gr2c,klev) 02527 case(126) 02528 !--- non-dimensionalize convective relaxation time scale 02529 if (mypid == NROOT) then 02530 zgp(:,:) = zgp(:,:)/ww 02531 endif 02532 allocate(gtdampc(nhor,klev)) 02533 if (mypid == NROOT) then 02534 write(nud,*) 'Field gtdampc allocated' 02535 endif 02536 call mpscgp(zgp,gtdampc,klev) 02537 end select 02538 kread = 1 02539 return 02540 end subroutine read_vargp 02541 02542 ! =============== 02543 ! SUBROUTINE DIAG 02544 ! =============== 02545 02546 subroutine diag 02547 use pumamod 02548 if (noutput > 0 .and. mod(nstep,ndiag) == 0) then 02549 if (ncoeff > 0) call prisp 02550 call xsect 02551 endif 02552 call energy 02553 return 02554 end 02555 02556 ! ================ 02557 ! SUBROUTINE PRISP 02558 ! ================ 02559 02560 subroutine prisp 02561 use pumamod 02562 02563 character(30) :: title 02564 02565 scale = 100.0 02566 title = 'Vorticity [10-2]' 02567 do 100 jlev=1,NLEV 02568 if (ndil(jlev).ne.0) call wrspam(sz(1,jlev),jlev,title,scale) 02569 100 continue 02570 02571 title = 'Divergence [10-2]' 02572 do 200 jlev=1,NLEV 02573 if (ndil(jlev).ne.0) call wrspam(sd(1,jlev),jlev,title,scale) 02574 200 continue 02575 02576 scale = 1000.0 02577 title = 'Temperature [10-3]' 02578 do 300 jlev=1,NLEV 02579 if (ndil(jlev).ne.0) call wrspam(st(1,jlev),jlev,title,scale) 02580 300 continue 02581 02582 title = 'Pressure [10-3]' 02583 call wrspam(sp,0,title,scale) 02584 02585 return 02586 end 02587 02588 ! ==================== 02589 ! SUBROUTINE POWERSPEC 02590 ! ==================== 02591 02592 subroutine powerspec(pf,pspec) 02593 use pumamod 02594 real :: pf(2,NCSP) 02595 real :: pspec(NTP1) 02596 02597 do j = 1 , NTP1 02598 pspec(j) = 0.5 * (pf(1,j) * pf(1,j) + pf(2,j) * pf(2,j)) 02599 enddo 02600 02601 j = NTP1 + 1 02602 do m = 2 , NTP1 02603 do l = m , NTP1 02604 pspec(l) = pspec(l) + pf(1,j) * pf(1,j) + pf(2,j) * pf(2,j) 02605 j = j + 1 02606 enddo 02607 enddo 02608 return 02609 end 02610 02611 ! ===================== 02612 ! SUBROUTINE POWERPRINT 02613 ! ===================== 02614 02615 subroutine powerprint(text,pspec) 02616 use pumamod 02617 character(3) :: text 02618 real :: pspec(NTP1) 02619 02620 zmax = maxval(pspec(:)) 02621 if (zmax <= 1.0e-20) return 02622 zsca = 10 ** (4 - int(log10(zmax))) 02623 write(nud,1000) text,(int(pspec(j)*zsca),j=2,13) 02624 return 02625 1000 format('* Power(',a3,') ',i8,11i5,' *') 02626 end 02627 02628 02629 02630 02631 ! ============== 02632 ! FUNCTION RMSSP 02633 ! ============== 02634 02635 function rmssp(pf) 02636 use pumamod 02637 real pf(NESP,NLEV) 02638 02639 zsum = 0.0 02640 do jlev = 1 , NLEV 02641 zsum = zsum + dsigma(jlev)& 02642 & * (dot_product(pf(1:NZOM,jlev),pf(1:NZOM,jlev)) * 0.5& 02643 & + dot_product(pf(NZOM+1:NRSP,jlev),pf(NZOM+1:NRSP,jlev))) 02644 enddo 02645 rmssp = zsum 02646 return 02647 end 02648 02649 ! ================= 02650 ! SUBROUTINE ENERGY 02651 ! ================= 02652 02653 subroutine energy 02654 use pumamod 02655 02656 parameter (idim=5) ! Number of scalars for GUI timeseries 02657 02658 ! calculates various global diagnostic quantities 02659 ! remove planetary vorticity so sz contains relative vorticity 02660 02661 real :: spec(NTP1) 02662 real (kind=4) ziso(idim) 02663 02664 sz(3,:) = sz(3,:) - plavor 02665 02666 ! *********************************************** 02667 ! calculate means - zpsitot rms vorticity 02668 ! zchitot rms divergence 02669 ! ztmptot rms temperature 02670 ! ztotp ie+pe potential energy 02671 ! zamsp mean surface pressure 02672 ! *********************************************** 02673 02674 zsqrt2 = sqrt(2.0) 02675 zamsp = 1.0 + span(1) / zsqrt2 02676 zst = dot_product(dsigma(:),st(1,:)) / zsqrt2 02677 ztout1 = dot_product(dsigma(:),t0(:)) 02678 02679 ztout2 = 0.0 02680 zst2b = 0.0 02681 ztoti = 0.0 02682 do jlev = 1 , NLEV 02683 ztout2 = ztout2 + dsigma(jlev) * t0(jlev) * t0(jlev) 02684 zst2b = zst2b + dsigma(jlev) * t0(jlev) * st(1,jlev) 02685 ztoti = ztoti + dsigma(jlev)& 02686 & * (dot_product(span(1:NZOM),st(1:NZOM,jlev)) * 0.5& 02687 & + dot_product(span(NZOM+1:NRSP),st(NZOM+1:NRSP,jlev))) 02688 enddo 02689 02690 ztotp = dot_product(span(1:NZOM),so(1:NZOM)) * 0.5& 02691 & + dot_product(span(NZOM+1:NRSP),so(NZOM+1:NRSP))& 02692 & + so(1)/zsqrt2 + (zamsp*ztout1+ztoti+zst) / akap 02693 02694 zpsitot = sqrt(rmssp(sz)) 02695 zchitot = sqrt(rmssp(sd)) 02696 ztmptot = sqrt(rmssp(st)+ztout2+zst2b*zsqrt2) 02697 02698 ziso(1) = ct * (spnorm(1) * st(1,NLEV) + t0(NLEV)) - 273.16 ! T(NLEV) [C] 02699 ziso(2) = ww * zchitot * 1.0e6 02700 ziso(3) = ztmptot 02701 ziso(4) = ztotp 02702 ziso(5) = sz(3,2) 02703 call guiput("SCALAR" // char(0) ,ziso,idim,1,1) 02704 02705 ! restore sz to absolute vorticity 02706 02707 sz(3,:) = sz(3,:) + plavor 02708 02709 if (mod(nstep,ndiag) /= 0) return ! was called for GUI only 02710 write(nud,9001) 02711 write(nud,9002) nstep,zpsitot,zchitot,ztmptot,ztotp,zamsp 02712 write(nud,9002) 02713 write(nud,9011) (j,j=1,12) 02714 write(nud,9012) 02715 call powerspec(span,spec) 02716 call powerprint('Pre',spec) 02717 call powerspec(sz(1,NLEV),spec) 02718 call powerprint('Vor',spec) 02719 call powerspec(sd(1,NLEV),spec) 02720 call powerprint('Div',spec) 02721 call powerspec(st(1,NLEV),spec) 02722 call powerprint('Tem',spec) 02723 return 02724 9001 format(/, 02725 ' nstep rms z rms d rms t & & pe+ie msp') 02726 9002 format(i10,4x,4g12.5,g15.8) 02727 !9009 format('*',75(' '),' *') 02728 !9010 format('* Power(',a,') ',7e9.2,' *') 02729 9011 format('* Wavenumber ',i8,11i5,' *') 02730 9012 format('',78('*')) 02731 end 02732 02733 ! ================= 02734 ! SUBROUTINE NTOMIN 02735 ! ================= 02736 02737 subroutine ntomin(kstep,imin,ihou,iday,imon,iyea) 02738 use pumamod 02739 istep = kstep ! day [0-29] month [0-11] 02740 if (istep .lt. 0) istep = 0 ! min [0-59] hour [0-23] 02741 imin = mod(istep,ntspd) * 1440 / ntspd ! minutes of current day 02742 ihou = imin / 60 ! hours of current day 02743 imin = imin - ihou * 60 ! minutes of current hour 02744 iday = istep / ntspd ! days in this run 02745 imon = iday / 30 ! months in this run 02746 iday = iday - imon * 30 ! days of current month 02747 iyea = imon / 12 ! years in this run 02748 imon = imon - iyea * 12 ! month of current year 02749 iday = iday + 1 02750 imon = imon + 1 02751 iyea = iyea + 1 02752 return 02753 end 02754 02755 ! ================= 02756 ! SUBROUTINE NTODAT 02757 ! ================= 02758 02759 subroutine ntodat(istep,datch) 02760 character(18) :: datch 02761 character(3) :: mona(12) 02762 data mona /'Jan','Feb','Mar','Apr','May','Jun',& 02763 & 'Jul','Aug','Sep','Oct','Nov','Dec'/ 02764 call ntomin(istep,imin,ihou,iday,imon,iyea) 02765 write(datch,20030) iday,mona(imon),iyea,ihou,imin 02766 20030 format(i2,'-',a3,'-',i4.4,2x,i2,':',i2.2) 02767 end 02768 02769 02770 ! ================= 02771 ! SUBROUTINE WRSPAM 02772 ! ================= 02773 02774 subroutine wrspam(ps,klev,title,scale) 02775 use pumamod 02776 ! 02777 dimension ps(NRSP) 02778 character(30) :: title 02779 character(18) :: datch 02780 02781 ! cab(i)=real(scale*sqrt(ps(i+i-1)*ps(i+i-1)+ps(i+i)*ps(i+i))) 02782 02783 call ntodat(nstep,datch) 02784 write(nud,'(1x)') 02785 write(nud,20000) 02786 write(nud,20030) datch,title,klev 02787 write(nud,20000) 02788 write(nud,20020) (i,i=0,9) 02789 write(nud,20000) 02790 write(nud,20100) (cab(i),i=1,10) 02791 write(nud,20200) (cab(i),i=NTRU+2,NTRU+10) 02792 write(nud,20300) (cab(i),i=2*NTRU+2,2*NTRU+9) 02793 write(nud,20400) (cab(i),i=3*NTRU+1,3*NTRU+7) 02794 write(nud,20000) 02795 write(nud,'(1x)') 02796 02797 20000 format(78('*')) 02798 20020 format('* n * ',10i7,' *') 02799 20030 format('* * ',a18,2x,a30,' Level ',i2,11x,'*') 02800 20100 format('* 0 *',f8.2,9f7.2,' *') 02801 20200 format('* 1 *',8x,9f7.2,' *') 02802 20300 format('* 2 *',15x,8f7.2,' *') 02803 20400 format('* 3 *',22x,7f7.2,' *') 02804 contains 02805 function cab(i) 02806 cab = scale * sqrt(ps(i+i-1)*ps(i+i-1)+ps(i+i)*ps(i+i)) 02807 end function cab 02808 end 02809 02810 ! =============== 02811 ! SUBROUTINE WRZS 02812 ! =============== 02813 02814 subroutine wrzs(zs,title,scale) 02815 use pumamod 02816 ! 02817 dimension zs(NLAT,NLEV) 02818 character(30) :: title 02819 character(18) :: datch 02820 02821 ip = NLAT / 16 02822 ia = ip/2 02823 ib = ia + 7 * ip 02824 id = NLAT + 1 - ia 02825 ic = id - 7 * ip 02826 02827 call ntodat(nstep,datch) 02828 write(nud,'(1x)') 02829 write(nud,20000) 02830 write(nud,20030) datch,title 02831 write(nud,20000) 02832 write(nud,20020) (chlat(i),i=ia,ib,ip),(chlat(j),j=ic,id,ip) 02833 write(nud,20000) 02834 do 200 jlev = 1 , NLEV 02835 write(nud,20100) jlev,((int(zs(i,jlev)*scale)),i=ia,ib,ip),& 02836 & ((int(zs(j,jlev)*scale)),j=ic,id,ip),jlev 02837 200 continue 02838 write(nud,20000) 02839 write(nud,'(1x)') 02840 02841 20000 format(78('*')) 02842 20020 format('* Lv * ',16(1x,a3),' * Lv *') 02843 20030 format('* * ',a18,2x,a30,20x,'*') 02844 20100 format('* ',i2,' * ',16i4,' * ',i2,' *') 02845 end 02846 02847 ! ================ 02848 ! SUBROUTINE XSECT 02849 ! ================ 02850 02851 subroutine xsect 02852 use pumamod 02853 character(30) :: title 02854 02855 scale = 10.0 02856 title = 'Zonal Wind [0.1 m/s]' 02857 call wrzs(csu,title,scale) 02858 title = 'Meridional Wind [0.1 m/s]' 02859 call wrzs(csv,title,scale) 02860 scale = 1.0 02861 title = 'Temperature [C]' 02862 call wrzs(cst,title,scale) 02863 return 02864 end 02865 02866 ! ================== 02867 ! SUBROUTINE WRITESP 02868 ! ================== 02869 02870 subroutine writesp(kunit,pf,kcode,klev,pscale,poff) 02871 use pumamod 02872 real :: pf(NRSP) 02873 real :: zf(NRSP) 02874 integer :: ihead(8) 02875 02876 call ntomin(nstep,nmin,nhour,nday,nmonth,nyear) 02877 02878 ihead(1) = kcode 02879 ihead(2) = klev 02880 ihead(3) = nday + 100 * nmonth + 10000 * nyear 02881 ihead(4) = nmin + 100 * nhour 02882 ihead(5) = NRSP 02883 ihead(6) = 1 02884 ihead(7) = 1 02885 ihead(8) = 0 02886 02887 ! normalize ECHAM compatible and scale to physical dimensions 02888 02889 zf(:) = pf(:) * spnorm(1:NRSP) * pscale 02890 zf(1) = zf(1) + poff ! Add offset if necessary 02891 write(kunit) ihead 02892 write(kunit) zf 02893 02894 return 02895 end 02896 02897 ! ================== 02898 ! SUBROUTINE WRITEGP 02899 ! ================== 02900 02901 subroutine writegp(kunit,pf,kcode,klev) 02902 use pumamod 02903 real :: pf(NHOR) 02904 real :: zf(NUGP) 02905 integer :: ihead(8) 02906 02907 call mpgagp(zf,pf,1) 02908 02909 if (mypid == NROOT) then 02910 call alt2reg(zf,1) 02911 call ntomin(nstep,nmin,nhour,nday,nmonth,nyear) 02912 02913 ihead(1) = kcode 02914 ihead(2) = klev 02915 ihead(3) = nday + 100 * nmonth + 10000 * nyear 02916 ihead(4) = nmin + 100 * nhour 02917 ihead(5) = NLON 02918 ihead(6) = NLAT 02919 ihead(7) = 1 02920 ihead(8) = 0 02921 02922 write(kunit) ihead 02923 write(kunit) zf 02924 endif 02925 02926 return 02927 end 02928 02929 02930 ! ================ 02931 ! SUBROUTINE OUTSP 02932 ! ================ 02933 02934 subroutine outsp 02935 use pumamod 02936 real zsr(NESP) 02937 02938 if (nwrioro == 1) then 02939 call writesp(40,so,129,0,cv*cv,0.0) 02940 nwrioro = 0 02941 endif 02942 02943 if (nextout == 1) then 02944 call writesp(40,sp2,40,0,1.0,log(psmean)) 02945 call writesp(40,sp1,41,0,1.0,log(psmean)) 02946 do jlev = 1,NLEV 02947 call writesp(40,st2(1,jlev),42,jlev,ct,t0(jlev)*ct) 02948 enddo 02949 do jlev = 1,NLEV 02950 call writesp(40,st1(1,jlev),43,jlev,ct,t0(jlev)*ct) 02951 enddo 02952 endif 02953 02954 ! ************ 02955 ! * pressure * 02956 ! ************ 02957 02958 call writesp(40,sp,152,0,1.0,log(psmean)) 02959 02960 ! *************** 02961 ! * temperature * 02962 ! *************** 02963 02964 do jlev = 1 , NLEV 02965 call writesp(40,st(1,jlev),130,jlev,ct,t0(jlev)*ct) 02966 enddo 02967 02968 ! ******************** 02969 ! * res. temperature * 02970 ! ******************** 02971 02972 zampl = cos((real(nstep)-pac)*tac) 02973 do jlev = 1 , NLEV 02974 zsr(:)=sr1(:,jlev)+sr2(:,jlev)*zampl 02975 call writesp(40,zsr,154,jlev,ct,t0(jlev)*ct) 02976 enddo 02977 02978 ! ************** 02979 ! * divergence * 02980 ! ************** 02981 02982 do jlev = 1 , NLEV 02983 call writesp(40,sd(1,jlev),155,jlev,ww,0.0) 02984 enddo 02985 02986 ! ************* 02987 ! * vorticity * 02988 ! ************* 02989 02990 do jlev = 1 , NLEV 02991 zsave = sz(3,jlev) 02992 sz(3,jlev) = sz(3,jlev) - plavor 02993 call writesp(40,sz(1,jlev),138,jlev,ww,0.0) 02994 sz(3,jlev) = zsave 02995 enddo 02996 02997 return 02998 end 02999 03000 ! ================ 03001 ! SUBROUTINE OUTGP 03002 ! ================ 03003 03004 subroutine outgp 03005 use pumamod 03006 real zhelp(NHOR) 03007 ! 03008 ! energy diagnostics 03009 ! 03010 if(nenergy > 0) then 03011 do je=1,9 03012 jcode=300+je 03013 zhelp(:)=denergy(:,je) 03014 call writegp(40,zhelp,jcode,0) 03015 enddo 03016 endif 03017 if(nentropy > 0) then 03018 do je=1,9 03019 jcode=310+je 03020 zhelp(:)=dentropy(:,je) 03021 call writegp(40,zhelp,jcode,0) 03022 enddo 03023 endif 03024 ! 03025 return 03026 end 03027 03028 03029 ! ==================== 03030 ! SUBROUTINE CHECKUNIT 03031 ! ==================== 03032 03033 subroutine checkunit 03034 use pumamod 03035 03036 write(ncu,1000) nstep,'sp( 1 )',sp(1),sp(1)*spnorm(1)+log(psmean) 03037 write(ncu,1000) nstep,'st( 1,1)',st(1,1),st(1,1)*spnorm(1)*ct+t0(1)*ct 03038 write(ncu,1000) nstep,'sd( 1,1)',sd(1,1),sd(1,1)*spnorm(1)*ww 03039 write(ncu,1000) nstep,'sz( 1,1)',sz(1,1),sz(1,1)*spnorm(1)*ww 03040 03041 write(ncu,1000) nstep,'st( 1,NLEV)',st(1,NLEV),st(1,NLEV)*spnorm(1)*ct+t0(5)*ct 03042 write(ncu,1000) nstep,'sd( 1,NLEV)',sd(1,NLEV),sd(1,NLEV)*spnorm(1)*ww 03043 write(ncu,1000) nstep,'sz( 1,NLEV)',sz(1,NLEV),sz(1,NLEV)*spnorm(1)*ww 03044 03045 if (100 < NRSP) then 03046 write(ncu,1000) nstep,'sp(100 )',sp(100),sp(100)*spnorm(100) 03047 write(ncu,1000) nstep,'st(100,NLEV)',st(100,NLEV),st(100,NLEV)*spnorm(100)*ct 03048 write(ncu,1000) nstep,'sd(100,NLEV)',sd(100,NLEV),sd(100,NLEV)*spnorm(100)*ww 03049 write(ncu,1000) nstep,'sz(100,NLEV)',sz(100,NLEV),sz(100,NLEV)*spnorm(100)*ww 03050 endif 03051 03052 return 03053 1000 format(i5,1x,a,1x,2f14.7) 03054 end 03055 03056 03057 ! ===================== 03058 ! * SUBROUTINE LEGPRI * 03059 ! ===================== 03060 03061 subroutine legpri 03062 use pumamod 03063 03064 write(nud,231) 03065 write(nud,232) 03066 write(nud,233) 03067 write(nud,232) 03068 do 14 jlat = 1 , NLAT 03069 zalat = asin(sid(jlat))*180.0/PI 03070 write(nud,234) jlat,zalat,csq(jlat),gwd(jlat) 03071 14 continue 03072 write(nud,232) 03073 write(nud,231) 03074 return 03075 231 format(/) 03076 232 format(37('*')) 03077 233 format('* No * Lat * csq weight *') 03078 234 format('*',i4,' *',f6.1,' *',2f10.4,' *') 03079 end 03080 03081 03082 ! ================= 03083 ! SUBROUTINE INILAT 03084 ! ================= 03085 03086 subroutine inilat 03087 use pumamod 03088 real (kind=8) :: zcsq 03089 03090 do jlat = 1 , NLAT 03091 zcsq = 1.0 - sid(jlat) * sid(jlat) 03092 csq(jlat) = zcsq 03093 rcs(jlat) = 1.0 / sqrt(zcsq) 03094 enddo 03095 do jlat = 1 , NLAT/2 03096 ideg = nint(180.0/PI * asin(sid(jlat))) 03097 write(chlat(jlat),'(i2,a1)') ideg,'N' 03098 write(chlat(NLAT+1-jlat),'(i2,a1)') ideg,'S' 03099 enddo 03100 return 03101 end 03102 03103 03104 ! ==================== 03105 ! SUBROUTINE GRIDPOINT 03106 ! ==================== 03107 03108 subroutine gridpoint 03109 use pumamod 03110 03111 real gtn(NLON,NLPP,NLEV) 03112 real gvpp(NHOR) 03113 real gpmt(NLON,NLPP) 03114 real sdf(NESP,NLEV) 03115 real stf(NESP,NLEV) 03116 real szf(NESP,NLEV) 03117 real spf(NESP) 03118 real zgp(NLON,NLAT) 03119 real zgpp(NHOR) 03120 real (kind=4) :: zcs(NLAT,NLEV) 03121 real (kind=4) :: zsp(NRSP) 03122 03123 do jlev = 1 , NLEV 03124 call sp2fc(sd(1,jlev),gd(1,jlev)) 03125 call sp2fc(st(1,jlev),gt(1,jlev)) 03126 call sp2fc(sz(1,jlev),gz(1,jlev)) 03127 enddo 03128 03129 call sp2fc(sp,gp) ! LnPs 03130 call sp2fcdmu(sp,gpj) ! d(lnps) / d(mu) 03131 ! divergence, vorticity -> u*cos(phi), v*cos(phi) 03132 do jlev = 1 , NLEV 03133 call dv2uv(sd(1,jlev),sz(1,jlev),gu(1,jlev),gv(1,jlev)) 03134 enddo 03135 if (lselect) then 03136 call filter_zonal_waves(gp) 03137 call filter_zonal_waves(gpj) 03138 do jlev = 1 , NLEV 03139 call filter_zonal_waves(gu(1,jlev)) 03140 call filter_zonal_waves(gv(1,jlev)) 03141 call filter_zonal_waves(gd(1,jlev)) 03142 call filter_zonal_waves(gt(1,jlev)) 03143 call filter_zonal_waves(gz(1,jlev)) 03144 enddo 03145 endif 03146 03147 if (ngui > 0 .or. mod(nstep,ndiag) == 0) then 03148 do jlev = 1 , NLEV 03149 do jlat = 1 , NLPP 03150 sec = cv / sqrt(csq(jlat)) 03151 csu(jlat,jlev) = gu(1+(jlat-1)*NLON,jlev) * sec 03152 csv(jlat,jlev) = gv(1+(jlat-1)*NLON,jlev) * sec 03153 cst(jlat,jlev) =(gt(1+(jlat-1)*NLON,jlev) + t0(jlev))*ct-273.16 03154 enddo 03155 enddo 03156 endif 03157 03158 do jlat = 1 , NLPP 03159 do jlon = 1 , NLON-1 , 2 03160 gpmt(jlon ,jlat) = -gp(jlon+1+(jlat-1)*NLON) * ((jlon-1)/2) 03161 gpmt(jlon+1,jlat) = gp(jlon +(jlat-1)*NLON) * ((jlon-1)/2) 03162 end do 03163 end do 03164 03165 call fc2gp(gu ,NLON,NLPP*NLEV) 03166 call fc2gp(gv ,NLON,NLPP*NLEV) 03167 call fc2gp(gt ,NLON,NLPP*NLEV) 03168 call fc2gp(gd ,NLON,NLPP*NLEV) 03169 call fc2gp(gz ,NLON,NLPP*NLEV) 03170 call fc2gp(gpj,NLON,NLPP) 03171 call fc2gp(gpmt,NLON,NLPP) 03172 03173 call calcgp(gtn,gpmt,gvpp) 03174 03175 gut(:,:) = gu(:,:) * gt(:,:) 03176 gvt(:,:) = gv(:,:) * gt(:,:) 03177 gke(:,:) = gu(:,:) * gu(:,:) + gv(:,:) * gv(:,:) 03178 03179 call gp2fc(gtn ,NLON,NLPP*NLEV) 03180 call gp2fc(gut ,NLON,NLPP*NLEV) 03181 call gp2fc(gvt ,NLON,NLPP*NLEV) 03182 call gp2fc(gfv ,NLON,NLPP*NLEV) 03183 call gp2fc(gfu ,NLON,NLPP*NLEV) 03184 call gp2fc(gke ,NLON,NLPP*NLEV) 03185 call gp2fc(gvpp,NLON,NLPP ) 03186 03187 call fc2sp(gvpp,spf) 03188 03189 if (lselect) then 03190 call filter_zonal_waves(gvpp) 03191 do jlev = 1 , NLEV 03192 call filter_zonal_waves(gtn(1,1,jlev)) 03193 call filter_zonal_waves(gut(1,jlev)) 03194 call filter_zonal_waves(gvt(1,jlev)) 03195 call filter_zonal_waves(gfv(1,jlev)) 03196 call filter_zonal_waves(gfu(1,jlev)) 03197 call filter_zonal_waves(gke(1,jlev)) 03198 enddo 03199 endif 03200 03201 do jlev = 1 , NLEV 03202 call mktend(sdf(1,jlev),stf(1,jlev),szf(1,jlev),gtn(1,1,jlev),& 03203 gfu(1,jlev),gfv(1,jlev),gke(1,jlev),gut(1,jlev),gvt(1,jlev)) 03204 enddo 03205 03206 if (nruido > 0) call stepruido 03207 call mpsumsc(spf,spt,1) 03208 call mpsumsc(stf,stt,NLEV) 03209 call mpsumsc(sdf,sdt,NLEV) 03210 call mpsumsc(szf,szt,NLEV) 03211 03212 if (ngui > 0 .or. mod(nstep,ndiag) == 0) then 03213 call fc2gp(gp,NLON,NLPP) 03214 zgpp(:) = exp(gp) ! LnPs -> Ps 03215 call mpgagp(zgp,zgpp,1) ! zgp = Ps (full grid) 03216 if (ngui > 0) then 03217 call guips(zgp,psmean) 03218 call guigv("GU" // char(0),gu) 03219 call guigv("GV" // char(0),gv) 03220 call guigt(gt) 03221 endif 03222 zgpp(:) = zgpp(:) - 1.0 ! Mean(LnPs) = 0 <-> Mean(Ps) = 1 03223 call gp2fc(zgpp,NLON,NLPP) 03224 call fc2sp(zgpp,span) 03225 03226 call mpsum(span,1) ! span = Ps spectral 03227 call mpgacs(csu) 03228 call mpgacs(csv) 03229 call mpgacs(cst) 03230 if (mypid == NROOT) then 03231 call altcs(csu) 03232 call altcs(csv) 03233 call altcs(cst) 03234 if (ngui > 0) then 03235 zcs(:,:) = csu(:,:) 03236 call guiput("CSU" // char(0) ,zcs ,NLAT,NLEV,1) 03237 zcs(:,:) = csv(:,:) 03238 call guiput("CSV" // char(0) ,zcs ,NLAT,NLEV,1) 03239 zcs(:,:) = cst(:,:) 03240 call guiput("CST" // char(0) ,zcs ,NLAT,NLEV,1) 03241 zsp(:) = span(1:NRSP) 03242 call guiput("SPAN" // char(0) ,zsp ,NCSP,-NTP1,1) 03243 endif 03244 endif 03245 endif 03246 return 03247 end 03248 03249 ! ================= 03250 ! SUBROUTINE CALCGP 03251 ! ================= 03252 subroutine calcgp(gtn,gpm,gvp) 03253 03254 use pumamod 03255 03256 ! Comments by Torben Kunz and Guido Schroeder 03257 03258 ! Compute nonlinear tendencies in grid point space. 03259 ! Hoskins and Simmons 1975 (Q.J.R.Meteorol.Soc.,101,637-655) (HS75) 03260 03261 ! For terms calculated in this routine, see HS75, eqs. (8)-(10) and 03262 ! appendix I: 03263 ! - script Fu, Fv as contributions to script D: gl. arrays gfu, gfv 03264 ! - script T: returned as gtn 03265 ! - script P: returned as gvp 03266 03267 03268 ! parameters (in) 03269 ! --------------- 03270 03271 ! gpm -- d(ln(ps)) / d(lambda) 03272 03273 ! parameters (out) 03274 ! --------------- 03275 03276 ! gtn -- temperature tendency 03277 ! gvp -- vertical integral of (u,v) * grad(ln(ps)) 03278 03279 ! global arrays variable in time 03280 ! ------------------------------ 03281 03282 ! gfu, gfv -- terms Fu, Fv in primitive equations, 03283 ! see HS75 (eqs. (1), (2)) 03284 ! gu, gv -- components u, v of horizontal velocity vector 03285 ! gd -- divergence D 03286 ! gz -- absolute vorticity 03287 ! gt -- temperature deviation T' 03288 03289 ! global arrays constant in time 03290 ! ------------------------------ 03291 03292 ! t0d -- reference temperature difference between two adjacent 03293 ! full levels 03294 ! tkp -- reference temperature times kappa (global parameter AKAP) 03295 ! rdsig -- 1 / (2 * dsigma) 03296 ! rcsq -- 1 / (1 - mu^2) 03297 03298 ! notations used in subsequent comments 03299 ! ------------------------------------- 03300 03301 ! aINTb(A)dsigma :<=> the integral of A over the interval [a,b] 03302 ! with respect to sigma 03303 03304 real gtn(NHOR,NLEV) 03305 real gpm(NHOR) , gvp(NHOR) 03306 real zsdotp(NHOR,NLEM),zsumd(NHOR),zsumvp(NHOR),zsumvpm(NHOR) 03307 real ztpta(NHOR),ztptb(NHOR) 03308 real zvgpg(NHOR,NLEV) 03309 real gtd(NHOR,NLEM) 03310 real gud(NHOR,NLEM) 03311 real gvd(NHOR,NLEM) 03312 03313 ! 1. 03314 ! 1.1 zvgpg: (u,v) * grad(ln(ps)) 03315 03316 do jlev = 1 , NLEV 03317 zvgpg(:,jlev) = rcsq * (gu(:,jlev)*gpm(:)+gv(:,jlev)*gpj(:)) 03318 enddo 03319 03320 ! 1.2 Calculate vertical integral of A = D + (u,v) * grad(ln(ps)), 03321 ! separated into divergence and ln(ps) advection. 03322 ! zsumd : 0INT1(D)dsigma 03323 ! gvp : 0INT1[(u,v) * grad ln(ps)]dsigma 03324 ! zsdotp : 0INTsigma(A)dsigma 03325 03326 zsumd = dsigma(1) * gd(:,1) 03327 gvp = dsigma(1) * zvgpg(:,1) 03328 zsdotp(:,1) = zsumd + gvp 03329 03330 do jlev = 2 , NLEM 03331 zsumd = zsumd + dsigma(jlev) * gd(:,jlev) 03332 gvp = gvp + dsigma(jlev) * zvgpg(:,jlev) 03333 zsdotp(:,jlev) = zsumd + gvp 03334 enddo 03335 03336 zsumd = zsumd + dsigma(NLEV) * gd(:,NLEV) 03337 gvp = gvp + dsigma(NLEV) * zvgpg(:,NLEV) 03338 03339 ! 2. Calculate vertical velocity and vertical advection terms 03340 ! on half levels. 03341 03342 do jlev = 1 , NLEM 03343 zsdotp(:,jlev) = (sigmh(jlev) * (zsumd+gvp) - zsdotp(:,jlev)) 03344 enddo 03345 03346 gtd(:,:) = zsdotp(:,:) * (gt(:,2:NLEV) - gt(:,1:NLEM)) 03347 gud(:,:) = zsdotp(:,:) * (gu(:,2:NLEV) - gu(:,1:NLEM)) 03348 gvd(:,:) = zsdotp(:,:) * (gv(:,2:NLEV) - gv(:,1:NLEM)) 03349 03350 ! 3. Calculate nonlinear contributions to temperature tendency and 03351 ! nonlinear terms Fu, Fv as used in vorticity and 03352 ! divergence equation. 03353 03354 ! 3.1 top level: 03355 03356 ! 3.1.1 zsumvp: 0INTsigma[(u,v) * grad(ln(ps))]dsigma 03357 03358 zsumvp = zvgpg(:,1) * dsigma(1) 03359 03360 ! 3.1.2 Calculation of gtn, gfv and gfu as for inner levels (3.2), 03361 ! but somewhat simplified: 03362 ! a) For the top level the following equation holds in the 03363 ! discretized form: (1/sigma)*0INTsigma(A)dsigma == A 03364 ! (HS75, second equation following eq. (7)). Therefore, 03365 ! (3.2.3) simplifies to -kappa*T' * D and (3.2.4) vanishes. 03366 ! b) Vertical advection terms (gtd, gud, gvd (see section 2) 03367 ! and vertical T0 advection (3.2.6)) vanish at upper 03368 ! boundary (sigma == 0). 03369 03370 gtn(:,1) = (1.0-akap) * gt(:,1) * gd(:,1) - rdsig(1) * (gtd(:,1) & 03371 + t0d(1) * (sigmh(1)*gvp-zsumvp)) 03372 03373 gfv(:,1) = - gu(:,1)*gz(:,1) - gpj(:)*gt(:,1) - rdsig(1)*gvd(:,1) 03374 gfu(:,1) = gv(:,1)*gz(:,1) - gpm(:)*gt(:,1) - rdsig(1)*gud(:,1) 03375 03376 ! 3.2 inner levels: 03377 03378 do jlev = 2 , NLEM 03379 03380 ! 3.2.1 ztpta: (1/sigma)*0INTsigma(A-D)dsigma 03381 ! ztptb: (1/sigma)*0INTsigma(A)dsigma 03382 ! Matrix c contains factors for discretized integration, see 03383 ! HS75 (second equation following eq. (7)). 03384 03385 ztpta = c(1,jlev) * zvgpg(:,1) 03386 ztptb = c(1,jlev) * (zvgpg(:,1) + gd(:,1)) 03387 03388 do jlej = 2 , jlev 03389 ztpta = ztpta + c(jlej,jlev) * zvgpg(:,jlej) 03390 ztptb = ztptb + c(jlej,jlev) * (zvgpg(:,jlej) + gd(:,jlej)) 03391 enddo 03392 03393 zsumvpm = zsumvp 03394 zsumvp = zsumvp + zvgpg(:,jlev) * dsigma(jlev) 03395 03396 ! 3.2.2 D * T' 03397 03398 gtn(:,jlev) = gt(:,jlev) * gd(:,jlev) 03399 03400 ! 3.2.3 kappa*T' * 03401 ! [(u,v)*grad(ln(ps)) - (1/sigma)*0INTsigma(A)dsigma] 03402 03403 gtn(:,jlev) = gtn(:,jlev) & 03404 & + akap * gt(:,jlev) * (zvgpg(:,jlev) - ztptb) 03405 03406 ! 3.2.4 kappa*T0 * 03407 ! [(u,v)*grad(ln(ps)) - (1/sigma)*0INTsigma(A-D)dsigma] 03408 03409 gtn(:,jlev) = gtn(:,jlev) & 03410 & + tkp(jlev) * (zvgpg(:,jlev) - ztpta) 03411 03412 ! 3.2.5 Calculate vertical T' advection on full levels by 03413 ! averaging two half level advection terms (gtd, calculated 03414 ! in section 2). 03415 03416 ! and 03417 03418 ! 3.2.6 Calculate vertical T0 advection on full levels by 03419 ! averaging two half level advection terms. 03420 03421 gtn(:,jlev) = gtn(:,jlev) & 03422 & - rdsig(jlev) * (gtd(:,jlev) + gtd(:,jlev-1) & 03423 & +(sigmh(jlev) * gvp - zsumvp) * t0d(jlev) & 03424 & +(sigmh(jlev-1) * gvp - zsumvpm) * t0d(jlev-1)) 03425 03426 ! 3.2.7 terms Fv, Fu, see HS75 (equations following eq. (5)); 03427 ! vertical advection terms interpolated to full levels by 03428 ! averaging two half level advection terms. 03429 03430 gfv(:,jlev) = - gu(:,jlev)*gz(:,jlev) - gpj(:)*gt(:,jlev) & 03431 & - rdsig(jlev)*(gvd(:,jlev) + gvd(:,jlev-1)) 03432 03433 gfu(:,jlev) = gv(:,jlev)*gz(:,jlev) - gpm(:)*gt(:,jlev) & 03434 & - rdsig(jlev)*(gud(:,jlev) + gud(:,jlev-1)) 03435 enddo 03436 03437 ! 3.3 bottom level 03438 03439 ! 3.3.1 ztpta, ztptb: see 3.2.1 03440 03441 ztpta = c(1,NLEV) * zvgpg(:,1) 03442 ztptb = c(1,NLEV) * (zvgpg(:,1) + gd(:,1)) 03443 03444 do jlej = 2 , NLEV 03445 ztpta = ztpta + c(jlej,NLEV) * zvgpg(:,jlej) 03446 ztptb = ztptb + c(jlej,NLEV) * (zvgpg(:,jlej) + gd(:,jlej)) 03447 enddo 03448 03449 ! 3.3.2 Calculation of gtn, gfv and gfu as for inner levels (3.2), 03450 ! but somewhat simplified: 03451 ! Vertical advection terms (gtd, gud, gvd (see section 2) and 03452 ! vertical T0 advection (3.2.6)) vanish at 03453 ! lower boundary (sigma == 1). 03454 03455 gtn(:,NLEV) = gt(:,NLEV) * gd(:,NLEV) & 03456 & + akap*gt(:,NLEV)*(zvgpg(:,NLEV)-ztptb) & 03457 & + tkp(NLEV)*(zvgpg(:,NLEV)-ztpta) & 03458 & - rdsig(NLEV) * (gtd(:,NLEM) & 03459 & + t0d(NLEM)*(sigmh(NLEM)*gvp-zsumvp)) 03460 03461 gfv(:,NLEV) = -gu(:,NLEV) * gz(:,NLEV) - gpj(:) * gt(:,NLEV) & 03462 & - rdsig(NLEV) * gvd(:,NLEM) 03463 gfu(:,NLEV) = gv(:,NLEV) * gz(:,NLEV) - gpm(:) * gt(:,NLEV) & 03464 & - rdsig(NLEV) * gud(:,NLEM) 03465 03466 ! 3.3.3 Add gaussian noise to T (controlled by nruido) 03467 03468 if (nruido > 0) gtn(:,:) = gtn(:,:) + ruidop(:,:) 03469 03470 return 03471 end 03472 03473 ! =================== 03474 ! SUBROUTINE SPECTRAL 03475 ! =================== 03476 03477 subroutine spectral 03478 use pumamod 03479 03480 !* Add adiabatic and diabatic tendencies - perform leapfrog 03481 03482 ! The adiabatic tendencies are added using the semi implicit scheme 03483 ! Hoskins and Simmons 1975 (Q.J.R.Meteorol.Soc.,101,637-655) (HS75) 03484 ! To compare the code directly with HS75 the following notes might 03485 ! be helpful (in addition to the comments below): 03486 03487 ! Name rule for global arrays <abc>: 03488 ! a : representation (s=spectral, g=grid, z=local) 03489 ! b : variable (p=ln(ps), d=divergence, z=vorticity, t=temperature) 03490 ! c : modifier (m=previous timestep, p=present timestep, t=tendency) 03491 03492 ! global arrays variable in time 03493 ! ------------------------------ 03494 03495 ! spt - pressure tendency HS75 (10) 03496 ! sdt - divergence tendency HS75 ( 8) 03497 ! szt - vorticity tendency 03498 ! stt - temperature tendency HS75 ( 9) 03499 03500 ! spm - pressure at previous timestep 03501 ! sdm - divergence at previous timestep 03502 ! szm - vorticity at previous timestep 03503 ! stm - temperature at previous timestep 03504 03505 ! spp - pressure at present timestep 03506 ! sdp - divergence at present timestep 03507 ! szp - vorticity at present timestep 03508 ! stp - temperature at present timestep 03509 03510 ! global arrays constant in time 03511 ! ------------------------------ 03512 03513 ! sak(NSPP) - = hyper diffusion 03514 ! sop(NSPP) - g* = orography as geopotential 03515 ! srp1(NSPP,NLEV) - Tr = radiative equilibrium temperature (annual mean) 03516 ! srp2(NSPP,NLEV) - Tr = radiative equilibrium temperature (annual cycle) 03517 ! nindex(NSPP) - n = total wavenumber n for spectral modes 03518 ! srcn(NSPP) - 1/Cn = 1.0 / (n * (n+1)) 03519 ! damp(NLEV) 1/tau R = time constant for newtonian cooling 03520 ! fric(NLEV) 1/tau F = time constant for Rayleigh friction 03521 03522 real zpm(NSPP) ! new spm 03523 real zdm(NSPP,NLEV) ! new sdm 03524 real zzm(NSPP,NLEV) ! new szm 03525 real ztm(NSPP,NLEV) ! new stm 03526 real zwp(NSPP) ! timefilter delta pm 03527 real zwd(NSPP,NLEV) ! timefilter delta sd 03528 real zwz(NSPP,NLEV) ! timefilter delta sz 03529 real zwt(NSPP,NLEV) ! timefilter delta st 03530 real zsrp(NSPP) ! restoring temperature (mean + annual cycle) 03531 03532 real zgt(NSPP,NLEV) ! work array 03533 03534 real,allocatable :: zstte(:,:,:) ! temp. tendencies for energy diag. 03535 real,allocatable :: zszte(:,:,:) ! vort. tendencies for energy recycling 03536 real,allocatable :: zsdte(:,:,:) ! div. tendencies for energy recycling 03537 real,allocatable :: zdps(:) ! surf pressure for energy diag 03538 real,allocatable :: zsp(:) ! surf pressure spectral 03539 real,allocatable :: zspf(:) ! surf pressure spectral 03540 real,allocatable :: zspt(:) ! surf pressure tendency 03541 03542 real,allocatable :: zst(:,:) ! temperature for entropy diagnostics 03543 real,allocatable :: zstt(:,:) ! tem. tendencies for entropy diag. 03544 real,allocatable :: ztgp(:,:) ! 03545 real,allocatable :: zdtgp(:,:) ! 03546 real,allocatable :: zsum1(:) 03547 real,allocatable :: zgw(:) 03548 03549 ! 0. Special code for experiments with mode filtering 03550 03551 if (lspecsel) call filter_spectral_modes 03552 03553 ! 1. Initialize local arrays 03554 03555 zpm(:) = spp(:) 03556 zdm(:,:) = sdp(:,:) 03557 zzm(:,:) = szp(:,:) 03558 ztm(:,:) = stp(:,:) 03559 ! 03560 ! allocate diagnostic arrays if needed 03561 ! 03562 if(nenergy > 0 .or. nentropy > 0) then 03563 allocate(zstte(NSPP,NLEV,3)) 03564 endif 03565 if(ndheat > 0) then 03566 allocate(zszte(NSPP,NLEV,2)) 03567 allocate(zsdte(NSPP,NLEV,2)) 03568 endif 03569 ! 03570 ! allocate and compute surface pressure if needed 03571 ! 03572 if(nenergy > 0 .or. nentropy > 0 .or. ndheat > 0) then 03573 allocate(zspt(NSPP)) 03574 allocate(zsp(NSPP)) 03575 endif 03576 03577 ! 2. Calculate divergence on timelevel t (sdt) HS75 (17) 03578 ! which will replace the divergence tendency sdt 03579 ! (semi implicit scheme) 03580 03581 ! The vertical scheme has being changed to the ECMWF scheme 03582 ! (see e.g. Simmons and Burridge 1981, Mon.Wea.Rev.,109,758-766). 03583 ! in this scheme, matrix xlphi (g) differs from that in HS75. 03584 03585 ! z0 : reference temperature To 03586 ! zq : 1.0 / Cn 03587 ! zt : xlphi * script T - To * script P 03588 ! zm : xlphi * T + To * ln(Ps)(t-dt) 03589 03590 ! (note that phi is needed in HS75 (17) and, therefore, 03591 ! the surface geopotential phi* [sop] is added 03592 03593 do jlev=1,NLEV 03594 z0 = t0(jlev) 03595 do jsp=1,NSPP 03596 zq = srcn(jsp) ! 1.0 / (n * (n + 1)) 03597 zt = dot_product(xlphi(:,jlev),stt(jsp,:)) - z0 * spt(jsp) 03598 zm = dot_product(xlphi(:,jlev),stm(jsp,:)) + z0 * spm(jsp) 03599 za = sdt(jsp,jlev) * zq 03600 zb = sdm(jsp,jlev) * zq 03601 zgt(jsp,jlev) = zb + delt * (za + zm + sop(jsp) + zt * delt) 03602 enddo 03603 enddo 03604 03605 ! bm1 is the invers of matrix (1/cn I+B dt**2) (lhs HS75 (17)) 03606 03607 do jlev = 1 , NLEV 03608 do jsp = 1 , NSPP 03609 jn = nindex(jsp) ! total wavenumber n 03610 sdt(jsp,jlev) = dot_product(zgt(jsp,:),bm1(:,jlev,jn)) 03611 enddo 03612 enddo 03613 03614 ! 3. Calculate surface pressure tendency -ln(ps) HS75 (15) 03615 03616 do jlev = 1 , NLEV 03617 spt(:) = spt(:) + dsigma(jlev) * sdt(:,jlev) 03618 enddo 03619 03620 ! 4. Calculate temperature tendency HS75 (14) 03621 03622 do jlev = 1 , NLEV 03623 do jsp = 1 , NSPP 03624 stt(jsp,jlev)=stt(jsp,jlev)-dot_product(xlt(:,jlev),sdt(jsp,:)) 03625 enddo 03626 enddo 03627 03628 ! 5. Add tendencies 03629 03630 spp(:) = spm(:) - delt2 * spt(:) ! spt = -ln(ps) tendency 03631 sdp(:,:) = 2.0 * sdt(:,:) - sdm(:,:) ! sdt = sdm + delt * tend. 03632 szp(:,:) = delt2 * szt(:,:) + szm(:,:) ! vorticity 03633 stp(:,:) = delt2 * stt(:,:) + stm(:,:) ! temperature 03634 03635 if(nenergy > 0) then 03636 zspt(:)=-spt(:) 03637 call mkenerdiag(stm,stt,spm,zspt,denergy(:,1)) 03638 endif 03639 if(nentropy > 0) then 03640 call mkentrodiag(stm,stt,spm,dentropy(:,1)) 03641 endif 03642 03643 ! 6. Calculate newtonian cooling, friction and biharmonic diffusion 03644 ! (srp - stp) * damp = (Tr' -T') / tau R = newtonian cooling 03645 ! srp1 = annual mean component 03646 ! srp2 = annual cycle component 03647 ! sak = diffusion 03648 ! fric = friction 03649 ! zampl = annual cycle 03650 03651 zampl = cos((real(nstep)-pac)*tac) 03652 03653 if (nhelsua == 0 .or. nhelsua == 1) then 03654 do jlev=1,NLEV 03655 zsrp(:)=srp1(:,jlev)+srp2(:,jlev)*zampl 03656 sdt(:,jlev) = sdp(:,jlev) * (sak(1:NSPP) - fric(jlev)) 03657 szt(:,jlev) = szp(:,jlev) * (sak(1:NSPP) - fric(jlev)) 03658 stt(:,jlev) = (zsrp(:) - stp(:,jlev)) * damp(jlev) & 03659 & + stp(:,jlev) * sak(1:NSPP) 03660 if(nenergy > 0) then 03661 zstte(:,jlev,2)=(zsrp(:)-stp(:,jlev))*damp(jlev) 03662 zstte(:,jlev,3)=stp(:,jlev)*sak(1:NSPP) 03663 endif 03664 if(ndheat > 0) then 03665 zsdte(:,jlev,1) = -sdp(:,jlev) * fric(jlev) 03666 zszte(:,jlev,1) = -szp(:,jlev) * fric(jlev) 03667 zsdte(:,jlev,2) = sdp(:,jlev) * sak(1:NSPP) 03668 zszte(:,jlev,2) = szp(:,jlev) * sak(1:NSPP) 03669 endif 03670 enddo 03671 elseif (nhelsua == 2 .or. nhelsua == 3 .or. ndiagp > 0) then 03672 if (ndiagp == 0) then 03673 call heatgp(zampl) ! stt(:,:) = Newtonian cooling 03674 else 03675 call diagp(zampl) ! stt(:,:) = Newtonian cooling 03676 endif 03677 if(nenergy > 0) then 03678 zstte(:,:,2)=stt(:,:) 03679 endif 03680 do jlev=1,NLEV 03681 sdt(:,jlev) = sdp(:,jlev) * (sak(1:NSPP) - fric(jlev)) 03682 szt(:,jlev) = szp(:,jlev) * (sak(1:NSPP) - fric(jlev)) 03683 stt(:,jlev) = stt(:,jlev) + stp(:,jlev) * sak(1:NSPP) 03684 if(nenergy > 0) then 03685 zstte(:,jlev,3)=stp(:,jlev)*sak(1:NSPP) 03686 endif 03687 if(ndheat > 0) then 03688 zsdte(:,jlev,1) = -sdp(:,jlev) * fric(jlev) 03689 zszte(:,jlev,1) = -szp(:,jlev) * fric(jlev) 03690 zsdte(:,jlev,2) = sdp(:,jlev) * sak(1:NSPP) 03691 zszte(:,jlev,2) = szp(:,jlev) * sak(1:NSPP) 03692 endif 03693 enddo 03694 endif 03695 03696 ! Conserve ln(ps) by forcing mode(0,0) to zero 03697 ! Correct vorticity by canceling the friction and diffusion 03698 ! applied to planetary vorticity 03699 ! Only root node processes the first NSPP modes 03700 03701 if(nenergy > 0) then 03702 zspt(:)=0. 03703 call mkenerdiag(stp,zstte(:,:,2),spp,zspt,denergy(:,2)) 03704 call mkenerdiag(stp,zstte(:,:,3),spp,zspt,denergy(:,3)) 03705 endif 03706 if(nentropy > 0) then 03707 call mkentrodiag(stp,zstte(:,:,2),spp,dentropy(:,2)) 03708 call mkentrodiag(stp,zstte(:,:,3),spp,dentropy(:,3)) 03709 endif 03710 if(nenergy > 0 .or. nentropy > 0 .or. ndheat > 0) then 03711 zsp(:)=spp(:) 03712 zstte(:,:,1)=stt(:,:) 03713 endif 03714 03715 if (mypid == NROOT) then 03716 spp(1) = 0.0 03717 spp(2) = 0.0 03718 szt(3,:) = szt(3,:) + plavor * (fric(:) - sak(3)) 03719 if(ndheat > 0) then 03720 zszte(3,:,1) = zszte(3,:,1) + plavor * fric(:) 03721 zszte(3,:,2) = zszte(3,:,2) - plavor * sak(3) 03722 endif 03723 endif 03724 ! 03725 ! 6b) call for vertical diffusion 03726 ! 03727 03728 if(dvdiff > 0.) call vdiff(stp,szp,sdp,stt,szt,sdt) 03729 03730 ! 03731 ! recycle kin energy dissipation 03732 ! 03733 03734 if(ndheat > 0) then 03735 call mkdheat(zszte(:,:,1),zszte(:,:,2) & 03736 & ,zsdte(:,:,1),zsdte(:,:,2),zsp) 03737 endif 03738 03739 03740 if(nenergy > 0 .or. nentropy > 0) then 03741 zstte(:,:,1)=stt(:,:)-zstte(:,:,1) 03742 endif 03743 if(nenergy > 0) then 03744 call mkenerdiag(stp,zstte(:,:,1),zsp,zspt,denergy(:,4)) 03745 endif 03746 if(nentropy > 0) then 03747 zstte(:,:,1)=stt(:,:)-zstte(:,:,1) 03748 call mkentrodiag(stp,zstte(:,:,1),zsp,dentropy(:,4)) 03749 endif 03750 if(nenergy > 0 .or. nentropy > 0) then 03751 zstte(:,:,1)=0. 03752 zspt(:)=(spp(:)-zsp(:))/delt2 03753 endif 03754 if(nenergy > 0) then 03755 call mkenerdiag(stp,zstte(:,:,1),zsp,zspt,denergy(:,8)) 03756 endif 03757 if(nentropy > 0) then 03758 call mkentrodiag(stp,zstte(:,:,1),zsp,dentropy(:,8)) 03759 endif 03760 03761 ! 03762 ! diagnostics of efficiency 03763 ! 03764 03765 if(ndheat > 1) then 03766 zcp=gascon/akap 03767 allocate(zst(NESP,NLEV)) 03768 allocate(zstt(NESP,NLEV)) 03769 allocate(zspf(NESP)) 03770 allocate(ztgp(NHOR,NLEV)) 03771 allocate(zdtgp(NHOR,NLEV)) 03772 allocate(zdps(NHOR)) 03773 allocate(zsum1(4)) 03774 allocate(zgw(NHOR)) 03775 jhor=0 03776 do jlat=1,NHPP 03777 do jlon=1,NLON*2 03778 jhor=jhor+1 03779 zgw(jhor)=gwd(jlat) 03780 enddo 03781 enddo 03782 call mpgallsp(zst,stp,NLEV) 03783 call mpgallsp(zstt,stt,NLEV) 03784 call mpgallsp(zspf,zsp,1) 03785 do jlev = 1 , NLEV 03786 call sp2fc(zst(1,jlev),ztgp(1,jlev)) 03787 call sp2fc(zstt(1,jlev),zdtgp(1,jlev)) 03788 enddo 03789 call sp2fc(zspf,zdps) 03790 call fc2gp(ztgp,NLON,NLPP*NLEV) 03791 call fc2gp(zdtgp,NLON,NLPP*NLEV) 03792 call fc2gp(zdps,NLON,NLPP) 03793 zdps(:)=psurf*exp(zdps(:)) 03794 zsum1(:)=0. 03795 do jlev=1,NLEV 03796 ztgp(:,jlev)=ct*(ztgp(:,jlev)+t0(jlev)) 03797 zdtgp(:,jlev)=ct*ww*zdtgp(:,jlev) 03798 zsum1(1)=zsum1(1)+SUM(zdtgp(:,jlev)*zgw(:) & 03799 & *zcp*zdps(:)/ga*dsigma(jlev) & 03800 & ,mask=(zdtgp(:,jlev) >= 0.)) 03801 zsum1(2)=zsum1(2)+SUM(zdtgp(:,jlev)*zgw(:) & 03802 & *zcp*zdps(:)/ga*dsigma(jlev) & 03803 & ,mask=(zdtgp(:,jlev) < 0.)) 03804 zsum1(3)=zsum1(3)+SUM(zdtgp(:,jlev)/ztgp(:,jlev)*zgw(:) & 03805 & *zcp*zdps(:)/ga*dsigma(jlev) & 03806 & ,mask=(zdtgp(:,jlev) >= 0.)) 03807 zsum1(4)=zsum1(4)+SUM(zdtgp(:,jlev)/ztgp(:,jlev)*zgw(:) & 03808 & *zcp*zdps(:)/ga*dsigma(jlev) & 03809 & ,mask=(zdtgp(:,jlev) < 0.)) 03810 enddo 03811 zsum3=SUM(zgw(:)) 03812 call mpsumbcr(zsum1,4) 03813 call mpsumbcr(zsum3,1) 03814 zsum1(:)=zsum1(:)/zsum3 03815 if(mypid == NROOT) then 03816 ztp=zsum1(1)/zsum1(3) 03817 zztm=zsum1(2)/zsum1(4) 03818 write(9,*) zsum1(:),zsum1(1)/zsum1(3),zsum1(2)/zsum1(4) & 03819 & ,(ztp-zztm)/ztp 03820 endif 03821 deallocate(zst) 03822 deallocate(zstt) 03823 deallocate(zspf) 03824 deallocate(ztgp) 03825 deallocate(zdps) 03826 deallocate(zdtgp) 03827 deallocate(zsum1) 03828 deallocate(zgw) 03829 endif 03830 03831 ! 7. Add newtonian cooling, friction and diffusion tendencies 03832 03833 sdp(:,:) = sdp(:,:) + delt2 * sdt(:,:) 03834 szp(:,:) = szp(:,:) + delt2 * szt(:,:) 03835 stp(:,:) = stp(:,:) + delt2 * stt(:,:) 03836 03837 ! 11. Coupling for synchronization runs 03838 03839 if (mrnum == 2 .and. nsync > 0) then 03840 call mrdiff(stp,std,NESP,NLEV) 03841 call mrdiff(sdp,sdd,NESP,NLEV) 03842 call mrdiff(szp,szd,NESP,NLEV) 03843 call mrdiff(spp,spd,NESP, 1) 03844 stp(:,:) = stp(:,:) + syncstr * std(:,:) 03845 sdp(:,:) = sdp(:,:) + syncstr * sdd(:,:) 03846 szp(:,:) = szp(:,:) + syncstr * szd(:,:) 03847 spp(: ) = spp(: ) + syncstr * spd(: ) 03848 03849 endif 03850 03851 ! 8. Apply Robert Asselin time filter (not for short initial timesteps) 03852 ! d(t) = pnu * f(t-1) + pnu * f(t+1) - 2 * pnu * f(t) 03853 03854 if (nkits == 0) then 03855 zwp(:) = pnu * (spm(:) + spp(:) - 2.0 * zpm(:) ) 03856 zwd(:,:) = pnu * (sdm(:,:) + sdp(:,:) - 2.0 * zdm(:,:)) 03857 zwz(:,:) = pnu * (szm(:,:) + szp(:,:) - 2.0 * zzm(:,:)) 03858 zwt(:,:) = pnu * (stm(:,:) + stp(:,:) - 2.0 * ztm(:,:)) 03859 03860 ! Add Robert-Asselin-Williams filter value to f(t) 03861 03862 spm(:) = zpm(:) + alpha * zwp(:) 03863 sdm(:,:) = zdm(:,:) + alpha * zwd(:,:) 03864 szm(:,:) = zzm(:,:) + alpha * zwz(:,:) 03865 stm(:,:) = ztm(:,:) + alpha * zwt(:,:) 03866 03867 ! Add filter value to f(t+1) 03868 03869 spp(:) = spp(:) - (1.0 - alpha) * zwp(:) 03870 sdp(:,:) = sdp(:,:) - (1.0 - alpha) * zwd(:,:) 03871 szp(:,:) = szp(:,:) - (1.0 - alpha) * zwz(:,:) 03872 stp(:,:) = stp(:,:) - (1.0 - alpha) * zwt(:,:) 03873 endif 03874 03875 if (nenergy > 0 .or. nentropy > 0) then 03876 zstte(:,:,1)=(stm(:,:)-ztm(:,:))/delt2 03877 zspt(:)=(spm(:)-zpm(:))/delt2 03878 endif 03879 if(nenergy > 0) then 03880 call mkenerdiag(ztm,zstte(:,:,1),zpm,zspt,denergy(:,9)) 03881 endif 03882 if (nentropy > 0) then 03883 call mkentrodiag(ztm,zstte(:,:,1),zpm,dentropy(:,9)) 03884 endif 03885 03886 ! 9. Save spectral arrays for extended output 03887 03888 if (nextout == 1 .and. mypid == NROOT) then 03889 if (mod(nstep,nafter) == nafter - 2) then 03890 if (.not. allocated(st2)) allocate(st2(nesp,nlev)) 03891 st2(:,:) = st(:,:) 03892 if (.not. allocated(sp2)) allocate(sp2(nesp)) 03893 sp2(:) = sp(:) 03894 endif 03895 if (mod(nstep,nafter) == nafter - 1) then 03896 if (.not. allocated(st1)) allocate(st1(nesp,nlev)) 03897 st1(:,:) = st(:,:) 03898 if (.not. allocated(sp1)) allocate(sp1(nesp)) 03899 sp1(:) = sp(:) 03900 endif 03901 endif 03902 03903 ! 10. Gather spectral modes from all processes 03904 03905 call mpgallsp(sp,spp, 1) 03906 call mpgallsp(sd,sdp,NLEV) 03907 call mpgallsp(sz,szp,NLEV) 03908 call mpgallsp(st,stp,NLEV) 03909 03910 if(nenergy > 0 .or. nentropy > 0) then 03911 deallocate(zstte) 03912 endif 03913 if(ndheat > 0) then 03914 deallocate(zszte) 03915 deallocate(zsdte) 03916 endif 03917 if(nenergy > 0 .or. nentropy > 0 .or. ndheat > 0) then 03918 deallocate(zsp) 03919 deallocate(zspt) 03920 endif 03921 03922 return 03923 end 03924 03925 subroutine mrcheck(f) 03926 use pumamod 03927 real :: f(*) 03928 write (nud,'(/,i3,8f8.4)') 0,f(1:16:2) 03929 write (nud,'( i3,8f8.4)') 8,f(17:32:2) 03930 write (nud,'( i3,8f8.4)') 16,f(33:48:2) 03931 write (nud,'( i3,8f8.4)') 24,f(49:64:2) 03932 write (nud,'( i3,8f8.4)') 32,f(65:80:2) 03933 return 03934 end 03935 03936 03937 ! ================ 03938 ! SUBROUTINE DIAGP 03939 ! ================ 03940 03941 subroutine diagp(zampl) 03942 use pumamod 03943 03944 real :: zstf(NESP,NLEV) 03945 real :: zgr12(NHOR,NLEV) 03946 real :: zgtt(NHOR,NLEV) 03947 real :: gr12(NHOR,NLEV) 03948 real :: gr12c(NHOR,NLEV) 03949 03950 03951 real :: gdtmp(NHOR) 03952 03953 real :: zampl 03954 03955 !--- transform temperature and divergence to grid point space 03956 call mpgallsp(st,stp,NLEV) 03957 if (nconv > 0) then 03958 call mpgallsp(sd,sdp,NLEV) 03959 endif 03960 do jlev=1,NLEV 03961 call sp2fc(st(1,jlev) ,gt(1,jlev) ) 03962 if (nconv > 0) then 03963 call sp2fc(sd(1,jlev) ,gd(1,jlev) ) 03964 endif 03965 enddo 03966 call fc2gp(gt ,NLON,NLPP*NLEV) 03967 if (nconv > 0) then 03968 call fc2gp(gd ,NLON,NLPP*NLEV) 03969 endif 03970 03971 03972 !--- radiative temperature tendencies 03973 gr12(:,:) = gr1(:,:) + gr2(:,:)*zampl 03974 zgtt(:,:) = (gr12(:,:) - gt(:,:))*gtdamp(:,:) 03975 03976 !--- add convective temperature tendencies 03977 if (nconv > 0) then 03978 gdtmp(:) = gd(:,nlev) 03979 do jlev = 1,nlev 03980 where (gdtmp < 0.0) 03981 gr12c(:,jlev) = gr1c(:,jlev) + gr2c(:,jlev)*zampl 03982 zgtt(:,jlev) = zgtt(:,jlev) + (gr12c(:,jlev) - gt(:,jlev))*gtdampc(:,jlev) 03983 endwhere 03984 enddo 03985 endif 03986 03987 !--- transform temperature tendencies to spectral space 03988 call gp2fc(zgtt ,NLON,NLPP*NLEV) 03989 do jlev=1,NLEV 03990 call fc2sp(zgtt(1,jlev),zstf(1,jlev)) 03991 enddo 03992 call mpsumsc(zstf,stt,NLEV) 03993 03994 return 03995 end subroutine diagp 03996 03997 ! ================= 03998 ! SUBROUTINE HEATGP 03999 ! ================= 04000 04001 subroutine heatgp(zampl) 04002 use pumamod 04003 04004 real :: zsr12(NESP,NLEV) 04005 real :: zsrp12(NSPP,NLEV) 04006 real :: zstf(NESP,NLEV) 04007 real :: zgr12(NHOR,NLEV) 04008 real :: zgtt(NHOR,NLEV) 04009 04010 real :: zampl 04011 04012 zsrp12(:,:)=srp1(:,:)+srp2(:,:)*zampl 04013 call mpgallsp(zsr12,zsrp12,NLEV) 04014 call mpgallsp(st,stp,NLEV) 04015 do jlev=1,NLEV 04016 call sp2fc(zsr12(1,jlev),zgr12(1,jlev)) 04017 call sp2fc(st(1,jlev) ,gt(1,jlev) ) 04018 enddo 04019 call fc2gp(zgr12,NLON,NLPP*NLEV) 04020 call fc2gp(gt ,NLON,NLPP*NLEV) 04021 04022 ! Newtonian cooling 04023 04024 zgtt(:,:) = (zgr12(:,:) - gt(:,:)) * gtdamp(:,:) 04025 04026 call gp2fc(zgtt ,NLON,NLPP*NLEV) 04027 do jlev=1,NLEV 04028 call fc2sp(zgtt(1,jlev),zstf(1,jlev)) 04029 enddo 04030 call mpsumsc(zstf,stt,NLEV) 04031 04032 return 04033 end 04034 04035 ! ================ 04036 ! SUBROUTINE VDIFF 04037 ! ================ 04038 04039 subroutine vdiff(pt,pz,pd,ptt,pzt,pdt) 04040 use pumamod 04041 ! 04042 parameter(ztref=250.) 04043 real pt(NSPP,NLEV),pz(NSPP,NLEV),pd(NSPP,NLEV) 04044 real ptt(NSPP,NLEV),pzt(NSPP,NLEV),pdt(NSPP,NLEV) 04045 real ztn(NSPP,NLEV),zzn(NSPP,NLEV),zdn(NSPP,NLEV) 04046 real zebs(NLEM) 04047 real zskap(NLEV),zskaph(NLEV) 04048 real zkdiff(NLEM) 04049 ! 04050 zdelt=delt2/ww 04051 zkonst1=ga*zdelt/gascon 04052 zkonst2=zkonst1*ga/gascon 04053 ! 04054 zskap(:)=sigma(:)**akap 04055 zskaph(:)=sigmh(:)**akap 04056 ! 04057 ! 1) modified diffusion coefficents 04058 ! 04059 do jlev=1,NLEM 04060 jlp=jlev+1 04061 zkdiff(jlev)=zkonst2*sigmh(jlev)*sigmh(jlev)/(ztref*ztref) & 04062 & *dvdiff/(sigma(jlp)-sigma(jlev)) 04063 enddo 04064 ! 04065 ! 2. semi implicit scheme 04066 ! 04067 ! 2a momentum 04068 ! 04069 ! top layer elimination 04070 ! 04071 zebs(1)=zkdiff(1)/(dsigma(1)+zkdiff(1)) 04072 zdn(:,1)=dsigma(1)*pd(:,1)/(dsigma(1)+zkdiff(1)) 04073 zzn(:,1)=dsigma(1)*pz(:,1)/(dsigma(1)+zkdiff(1)) 04074 ! 04075 ! middle layer elimination 04076 ! 04077 do jlev=2,NLEM 04078 jlem=jlev-1 04079 zebs(jlev)=zkdiff(jlev)/(dsigma(jlev)+zkdiff(jlev) & 04080 & +zkdiff(jlem)*(1.-zebs(jlem))) 04081 zdn(:,jlev)=(pd(:,jlev)*dsigma(jlev)+zkdiff(jlem)*zdn(:,jlem)) & 04082 & /(dsigma(jlev)+zkdiff(jlev) & 04083 & +zkdiff(jlem)*(1.-zebs(jlem))) 04084 zzn(:,jlev)=(pz(:,jlev)*dsigma(jlev)+zkdiff(jlem)*zzn(:,jlem)) & 04085 & /(dsigma(jlev)+zkdiff(jlev) & 04086 & +zkdiff(jlem)*(1.-zebs(jlem))) 04087 enddo 04088 ! 04089 ! bottom layer elimination 04090 ! 04091 zdn(:,NLEV)=(pd(:,NLEV)*dsigma(NLEV)+zkdiff(NLEM)*zdn(:,NLEM)) & 04092 & /(dsigma(NLEV)+zkdiff(NLEM)*(1.-zebs(NLEM))) 04093 zzn(:,NLEV)=(pz(:,NLEV)*dsigma(NLEV)+zkdiff(NLEM)*zzn(:,NLEM)) & 04094 & /(dsigma(NLEV)+zkdiff(NLEM)*(1.-zebs(NLEM))) 04095 ! 04096 ! back-substitution 04097 ! 04098 do jlev=NLEM,1,-1 04099 jlep=jlev+1 04100 zdn(:,jlev)=zdn(:,jlev)+zebs(jlev)*zdn(:,jlep) 04101 zzn(:,jlev)=zzn(:,jlev)+zebs(jlev)*zzn(:,jlep) 04102 enddo 04103 ! 04104 ! tendencies 04105 ! 04106 pdt(:,1:NLEV)=pdt(:,1:NLEV)+(zdn(:,1:NLEV)-pd(:,1:NLEV))/delt2 04107 pzt(:,1:NLEV)=pzt(:,1:NLEV)+(zzn(:,1:NLEV)-pz(:,1:NLEV))/delt2 04108 ! 04109 ! 2c potential temperature 04110 ! 04111 do jlev=1,NLEM 04112 zkdiff(jlev)=zkdiff(jlev)*zskaph(jlev) 04113 enddo 04114 ! 04115 ! semi implicit scheme 04116 ! 04117 ! top layer elimination 04118 ! 04119 zebs(1)=zkdiff(1)/(dsigma(1)+zkdiff(1)/zskap(1)) 04120 ztn(:,1)=dsigma(1)*pt(:,1)/(dsigma(1)+zkdiff(1)/zskap(1)) 04121 ! 04122 ! middle layer elimination 04123 ! 04124 do jlev=2,NLEM 04125 jlem=jlev-1 04126 zebs(jlev)=zkdiff(jlev)/(dsigma(jlev)+(zkdiff(jlev) & 04127 & +zkdiff(jlem)*(1.-zebs(jlem)/zskap(jlem)))/zskap(jlev)) 04128 ztn(:,jlev)=(pt(:,jlev)*dsigma(jlev) & 04129 & +zkdiff(jlem)/zskap(jlem)*ztn(:,jlem)) & 04130 & /(dsigma(jlev)+(zkdiff(jlev) & 04131 & +zkdiff(jlem)*(1.-zebs(jlem)/zskap(jlem))) & 04132 & /zskap(jlev)) 04133 enddo 04134 ! 04135 ! bottom layer elimination 04136 ! 04137 ztn(:,NLEV)=(pt(:,NLEV)*dsigma(NLEV) & 04138 & +zkdiff(NLEM)*ztn(:,NLEM)/zskap(NLEM)) & 04139 & /(dsigma(NLEV)+zkdiff(NLEM)/zskap(NLEV) & 04140 & *(1.-zebs(NLEM)/zskap(NLEM))) 04141 ! 04142 ! back-substitution 04143 ! 04144 do jlev=NLEM,1,-1 04145 jlep=jlev+1 04146 ztn(:,jlev)=ztn(:,jlev)+zebs(jlev)*ztn(:,jlep)/zskap(jlep) 04147 enddo 04148 ! 04149 ! tendencies 04150 ! 04151 ptt(:,1:NLEV)=ptt(:,1:NLEV)+(ztn(:,1:NLEV)-pt(:,1:NLEV))/delt2 04152 ! 04153 return 04154 end subroutine vdiff 04155 04156 ! ================= 04157 ! SUBROUTINE GASDEV 04158 ! ================= 04159 04160 ! Gaussian noise generator with zero mean and unit variance. 04161 04162 real function gasdev() 04163 use pumamod 04164 implicit none 04165 real :: fr, vx, vy, ra 04166 04167 if (ganext == 0.0) then 04168 ra = 2.0 04169 do while (ra >= 1.0 .or. ra < 1.0e-20) 04170 call random_number(vx) 04171 call random_number(vy) 04172 vx = 2.0 * vx - 1.0 04173 vy = 2.0 * vy - 1.0 04174 ra = vx * vx + vy * vy 04175 enddo 04176 fr = sqrt(-2.0 * log(ra) / ra) 04177 gasdev = vx * fr 04178 ganext = vy * fr 04179 else 04180 gasdev = ganext 04181 ganext = 0.0 04182 endif 04183 04184 return 04185 end 04186 04187 ! ================= 04188 ! SUBROUTINE SPONGE 04189 ! ================= 04190 04191 subroutine sponge 04192 use pumamod 04193 04194 real :: zp 04195 04196 ! This introduces a simple sponge layer to the highest model levels 04197 ! by applying Rayleigh friction there, according to 04198 ! Polvani & Kushner (2002, GRL), see their appendix. 04199 04200 write(nud,*) 04201 write(nud,9991) 04202 write(nud,9997) 04203 write(nud,9991) 04204 write(nud,9996) 04205 write(nud,9991) 04206 do jlev=1,NLEV 04207 zp = sigma(jlev)*psurf 04208 if (zp < pspon) then 04209 fric(jlev) = (sponk * ((pspon - zp) / pspon)**2) / TWOPI 04210 endif 04211 04212 ! some output 04213 if (zp > pspon) then 04214 if (fric(jlev) == 0) then 04215 write(nud,9992) jlev 04216 else 04217 write(nud,9993) jlev, fric(jlev)*TWOPI 04218 endif 04219 else 04220 if (fric(jlev) == 0) then 04221 write(nud,9994) jlev 04222 else 04223 write(nud,9995) jlev, fric(jlev)*TWOPI 04224 endif 04225 endif 04226 enddo 04227 write(nud,9991) 04228 write(nud,*) 04229 return 04230 9991 format(33('*')) 04231 9992 format('*',i4,' * ',7('-'),' * *') 04232 9993 format('*',i4,' * ',f7.4,' * *') 04233 9994 format('*',i4,' * ',7('-'),' *',' SPONGE *') 04234 9995 format('*',i4,' * ',f7.4,' *',' SPONGE *') 04235 9996 format('* Lv * [1/day] * *') 04236 9997 format('* Rayleigh damping coefficients *') 04237 end 04238 04239 04240 ! ===================== 04241 ! SUBROUTINE MKENERDIAG 04242 ! ===================== 04243 04244 subroutine mkenerdiag(pst,pstt,psp,pspt,penergy) 04245 use pumamod 04246 ! 04247 real :: pst(NSPP,NLEV),pstt(NSPP,NLEV) 04248 real :: psp(NSPP),pspt(NSPP) 04249 real :: penergy(NHOR) 04250 ! 04251 real :: zsttf(NESP,NLEV),zstf(NESP,NLEV) 04252 real :: zsptf(NESP),zspf(NESP) 04253 real :: zgtt(NHOR,NLEV),zgt(NHOR,NLEV) 04254 real :: zgps(NHOR),zgpst(NHOR) 04255 real :: ztm(NHOR) 04256 ! 04257 zcp=gascon/akap 04258 zdelt=delt2/ww 04259 ! 04260 call mpgallsp(zsttf,pstt,NLEV) 04261 call mpgallsp(zstf,pst,NLEV) 04262 call mpgallsp(zsptf,pspt,1) 04263 call mpgallsp(zspf,psp,1) 04264 04265 do jlev=1,NLEV 04266 call sp2fc(zsttf(:,jlev),zgtt(:,jlev)) 04267 call sp2fc(zstf(:,jlev),zgt(:,jlev)) 04268 enddo 04269 call sp2fc(zsptf,zgpst) 04270 call sp2fc(zspf,zgps) 04271 call fc2gp(zgtt,NLON,NLPP*NLEV) 04272 call fc2gp(zgt,NLON,NLPP*NLEV) 04273 call fc2gp(zgps,NLON,NLPP) 04274 call fc2gp(zgpst,NLON,NLPP) 04275 zgpst(:)=psurf*(exp(zgps(:)+delt2*zgpst(:))-exp(zgps(:)))/zdelt 04276 zgps(:)=psurf*exp(zgps(:))+zdelt*zgpst(:) 04277 zgtt(:,:)=ct*ww*zgtt(:,:) 04278 do jlev=1,NLEV 04279 zgt(:,jlev)=ct*(zgt(:,jlev)+t0(jlev)) 04280 enddo 04281 ! 04282 ztm(:)=0. 04283 penergy(:)=0. 04284 do jlev=1,NLEV 04285 ztm(:)=ztm(:)+zgt(:,jlev)*dsigma(jlev) 04286 penergy(:)=penergy(:)+zgtt(:,jlev)*dsigma(jlev) 04287 enddo 04288 penergy(:)=ztm(:)*zcp*zgpst(:)/ga & 04289 & +penergy(:)*zcp*zgps(:)/ga 04290 ! 04291 return 04292 end 04293 04294 ! ====================== 04295 ! SUBROUTINE MKENTRODIAG 04296 ! ====================== 04297 04298 subroutine mkentrodiag(pst,pstt,psp,pentropy) 04299 use pumamod 04300 ! 04301 real :: pst(NSPP,NLEV),pstt(NSPP,NLEV) 04302 real :: psp(NSPP) 04303 real :: pentropy(NHOR) 04304 ! 04305 real :: zsttf(NESP,NLEV),zstf(NESP,NLEV) 04306 real :: zspf(NESP) 04307 real :: zgtt(NHOR,NLEV),zgt(NHOR,NLEV) 04308 real :: zgps(NHOR) 04309 ! 04310 zcp=gascon/akap 04311 ! 04312 call mpgallsp(zsttf,pstt,NLEV) 04313 call mpgallsp(zstf,pst,NLEV) 04314 call mpgallsp(zspf,psp,1) 04315 04316 do jlev=1,NLEV 04317 call sp2fc(zsttf(:,jlev),zgtt(:,jlev)) 04318 call sp2fc(zstf(:,jlev),zgt(:,jlev)) 04319 enddo 04320 call sp2fc(zspf,zgps) 04321 call fc2gp(zgtt,NLON,NLPP*NLEV) 04322 call fc2gp(zgt,NLON,NLPP*NLEV) 04323 call fc2gp(zgps,NLON,NLPP) 04324 zgps(:)=psurf*exp(zgps(:)) 04325 zgtt(:,:)=ct*ww*zgtt(:,:) 04326 do jlev=1,NLEV 04327 zgt(:,jlev)=ct*(zgt(:,jlev)+t0(jlev)) 04328 enddo 04329 ! 04330 pentropy(:)=0. 04331 do jlev=1,NLEV 04332 pentropy(:)=pentropy(:)+zgtt(:,jlev)*dsigma(jlev)/zgt(:,jlev) 04333 enddo 04334 pentropy(:)=pentropy(:)*zcp*zgps(:)/ga 04335 ! 04336 return 04337 end 04338 04339 ! ================== 04340 ! SUBROUTINE MKDHEAT 04341 ! ================== 04342 04343 subroutine mkdheat(zszt1,zszt2,zsdt1,zsdt2,zsp) 04344 use pumamod 04345 ! 04346 ! 'recycle' kin. energy loss by heating the environment 04347 ! 04348 ! zszt1/zsdt1 : vorticity/divergence tendency due to friction 04349 ! zszt2/zsdt2 : vorticity/divergence tendency fue to diffusion 04350 ! zp : surface pressure 04351 ! 04352 real zszt1(NSPP,NLEV),zszt2(NSPP,NLEV) 04353 real zsdt1(NSPP,NLEV),zsdt2(NSPP,NLEV) 04354 real zsp(NSPP) 04355 real zp(NHOR) 04356 ! 04357 real zsd(NESP,NLEV),zsz(NESP,NLEV) 04358 real zspf(NESP),zspt(NSPP) 04359 real zsdp(NSPP,NLEV),zszp(NSPP,NLEV) 04360 real zu(NHOR,NLEV),zun(NHOR,NLEV),zdu1(NHOR,NLEV),zdu2(NHOR,NLEV) 04361 real zv(NHOR,NLEV),zvn(NHOR,NLEV),zdv1(NHOR,NLEV),zdv2(NHOR,NLEV) 04362 real zdtdt1(NHOR,NLEV),zdtdt2(NHOR,NLEV),zdtdt3(NHOR,NLEV) 04363 ! 04364 real zdtdt(NHOR,NLEV),zdekin(NHOR,NLEV) 04365 ! 04366 real zsde(NSPP,NLEV),zsdef(NESP,NLEV) 04367 real zstt(NSPP,NLEV),zstf(NESP,NLEV) 04368 real zstt1(NSPP,NLEV),zstf1(NESP,NLEV),zstt3(NSPP,NLEV) 04369 real zstt2(NSPP,NLEV),zstf2(NESP,NLEV),zstf3(NESP,NLEV) 04370 ! 04371 ! some constants 04372 ! 04373 zdelt=delt2/ww ! timestep in s 04374 zcp=gascon/akap ! heat capacity 04375 ! 04376 ! 'recycle' friction 04377 ! 04378 ! a) gather the 'partial' field of z and d, and make u and v 04379 ! at old time level 04380 ! 04381 zsdp(:,:)=sdp(:,:) 04382 zszp(:,:)=szp(:,:) 04383 call mpgallsp(zsd,zsdp,NLEV) 04384 call mpgallsp(zsz,zszp,NLEV) 04385 do jlev = 1 , NLEV 04386 call dv2uv(zsd(1,jlev),zsz(1,jlev),zu(1,jlev),zv(1,jlev)) 04387 enddo 04388 call fc2gp(zu,NLON,NLPP*NLEV) 04389 call fc2gp(zv,NLON,NLPP*NLEV) 04390 ! 04391 ! b) add fricton tendencies and create new u and v 04392 ! 04393 zsdp(:,:)=sdp(:,:)+zsdt1(:,:)*delt2 04394 zszp(:,:)=szp(:,:)+zszt1(:,:)*delt2 04395 call mpgallsp(zsd,zsdp,NLEV) 04396 call mpgallsp(zsz,zszp,NLEV) 04397 do jlev = 1 , NLEV 04398 call dv2uv(zsd(1,jlev),zsz(1,jlev),zun(1,jlev),zvn(1,jlev)) 04399 enddo 04400 call fc2gp(zun,NLON,NLPP*NLEV) 04401 call fc2gp(zvn,NLON,NLPP*NLEV) 04402 ! 04403 ! c) compute temperature tendency 04404 ! 04405 do jlev=1,NLEV 04406 zu(:,jlev)=cv*zu(:,jlev)*SQRT(rcsq(:)) 04407 zv(:,jlev)=cv*zv(:,jlev)*SQRT(rcsq(:)) 04408 zun(:,jlev)=cv*zun(:,jlev)*SQRT(rcsq(:)) 04409 zvn(:,jlev)=cv*zvn(:,jlev)*SQRT(rcsq(:)) 04410 zdu1(:,jlev)=zun(:,jlev)-zu(:,jlev) 04411 zdv1(:,jlev)=zvn(:,jlev)-zv(:,jlev) 04412 zdtdt1(:,jlev)=-(zun(:,jlev)*zun(:,jlev) & 04413 & -zu(:,jlev)*zu(:,jlev) & 04414 & +zvn(:,jlev)*zvn(:,jlev) & 04415 & -zv(:,jlev)*zv(:,jlev))*0.5/zdelt/zcp 04416 enddo 04417 04418 ! 04419 ! 'recycle' momentum diffusion 04420 ! 04421 ! a) add tendencies and create new u and v and get surface pressure 04422 ! 04423 ! 04424 zsdp(:,:)=sdp(:,:)+zsdt2(:,:)*delt2 04425 zszp(:,:)=szp(:,:)+zszt2(:,:)*delt2 04426 call mpgallsp(zsd,zsdp,NLEV) 04427 call mpgallsp(zsz,zszp,NLEV) 04428 call mpgallsp(zspf,zsp,1) 04429 do jlev = 1 , NLEV 04430 call dv2uv(zsd(1,jlev),zsz(1,jlev),zun(1,jlev),zvn(1,jlev)) 04431 enddo 04432 call fc2gp(zun,NLON,NLPP*NLEV) 04433 call fc2gp(zvn,NLON,NLPP*NLEV) 04434 call sp2fc(zspf,zp) 04435 call fc2gp(zp,NLON,NLPP) 04436 zp(:)=psurf*exp(zp(:)) 04437 ! 04438 ! b) compute loss of kinetic energy 04439 ! (note: only the global average change of kin. e. is 'lost' 04440 ! the other changes are just diffusion) 04441 ! 04442 do jlev = 1 , NLEV 04443 zun(:,jlev)=cv*zun(:,jlev)*SQRT(rcsq(:)) 04444 zvn(:,jlev)=cv*zvn(:,jlev)*SQRT(rcsq(:)) 04445 zdu2(:,jlev)=zun(:,jlev)-zu(:,jlev) 04446 zdv2(:,jlev)=zvn(:,jlev)-zv(:,jlev) 04447 zdekin(:,jlev)=(zun(:,jlev)*zun(:,jlev) & 04448 & -zu(:,jlev)*zu(:,jlev) & 04449 & +zvn(:,jlev)*zvn(:,jlev) & 04450 & -zv(:,jlev)*zv(:,jlev))*0.5/zdelt & 04451 & *zp(:)/ga*dsigma(jlev) 04452 enddo 04453 ! 04454 ! c) get the global average and transform it back 04455 ! 04456 call gp2fc(zdekin,NLON,NLPP*NLEV) 04457 do jlev=1,NLEV 04458 call fc2sp(zdekin(:,jlev),zsdef(:,jlev)) 04459 enddo 04460 call mpsumsc(zsdef,zsde,NLEV) 04461 call mpgallsp(zsdef,zsde,NLEV) 04462 zsdef(2:NESP,:)=0. 04463 do jlev = 1 , NLEV 04464 call sp2fc(zsdef(1,jlev),zdekin(1,jlev)) 04465 enddo 04466 call fc2gp(zdekin,NLON,NLPP*NLEV) 04467 ! 04468 ! d) compute temperature tendency 04469 ! 04470 do jlev=1,NLEV 04471 zdtdt2(:,jlev)=-zdekin(:,jlev)*ga/zp(:)/dsigma(jlev)/zcp 04472 zdtdt3(:,jlev)=-(zdu1(:,jlev)*zdu2(:,jlev) & 04473 & +zdv1(:,jlev)*zdv2(:,jlev))/zdelt/zcp 04474 enddo 04475 ! 04476 zdtdt1(:,:)=zdtdt1(:,:)/ct/ww 04477 zdtdt2(:,:)=zdtdt2(:,:)/ct/ww 04478 zdtdt3(:,:)=zdtdt3(:,:)/ct/ww 04479 ! 04480 call gp2fc(zdtdt1,NLON,NLPP*NLEV) 04481 call gp2fc(zdtdt2,NLON,NLPP*NLEV) 04482 call gp2fc(zdtdt3,NLON,NLPP*NLEV) 04483 do jlev=1,NLEV 04484 call fc2sp(zdtdt1(:,jlev),zstf1(:,jlev)) 04485 call fc2sp(zdtdt2(:,jlev),zstf2(:,jlev)) 04486 call fc2sp(zdtdt3(:,jlev),zstf3(:,jlev)) 04487 enddo 04488 call mpsumsc(zstf1,zstt1,NLEV) 04489 call mpsumsc(zstf2,zstt2,NLEV) 04490 call mpsumsc(zstf3,zstt3,NLEV) 04491 ! 04492 ! add the temprature tendencies 04493 ! 04494 stt(:,:)=stt(:,:)+zstt1(:,:)+zstt2(:,:)+zstt3(:,:) 04495 ! 04496 ! energy diagnostics 04497 ! 04498 if(nenergy > 0) then 04499 zspt(:)=0. 04500 call mkenerdiag(stp,zstt1,zsp,zspt,denergy(:,5)) 04501 call mkenerdiag(stp,zstt2,zsp,zspt,denergy(:,6)) 04502 call mkenerdiag(stp,zstt3,zsp,zspt,denergy(:,7)) 04503 endif 04504 if(nentropy > 0) then 04505 call mkentrodiag(stp,zstt1,zsp,dentropy(:,5)) 04506 call mkentrodiag(stp,zstt2,zsp,dentropy(:,6)) 04507 call mkentrodiag(stp,zstt3,zsp,dentropy(:,7)) 04508 endif 04509 04510 ! 04511 return 04512 end subroutine mkdheat 04513 04514 ! ================= 04515 ! SUBROUTINE MKEKIN 04516 ! ================= 04517 04518 subroutine mkekin(zszp,zsdp,zp,zekin) 04519 use pumamod 04520 ! 04521 real zszp(NSPP,NLEV),zsdp(NSPP,NLEV) 04522 real zp(NHOR),zekin(NHOR) 04523 ! 04524 real zsd(NESP,NLEV),zsz(NESP,NLEV) 04525 real zu(NHOR,NLEV),zv(NHOR,NLEV) 04526 ! 04527 ! some constants 04528 ! 04529 zdelt=delt2/ww ! timestep in s 04530 zcp=gascon/akap ! heat capacity 04531 ! 04532 call mpgallsp(zsd,zsdp,NLEV) 04533 call mpgallsp(zsz,zszp,NLEV) 04534 do jlev = 1 , NLEV 04535 call dv2uv(zsd(1,jlev),zsz(1,jlev),zu(1,jlev),zv(1,jlev)) 04536 enddo 04537 call fc2gp(zu,NLON,NLPP*NLEV) 04538 call fc2gp(zv,NLON,NLPP*NLEV) 04539 ! 04540 zekin(:)=0. 04541 do jlev = 1 , NLEV 04542 zu(:,jlev)=cv*zu(:,jlev)*SQRT(rcsq(:)) 04543 zv(:,jlev)=cv*zv(:,jlev)*SQRT(rcsq(:)) 04544 zekin(:)=(zu(:,jlev)*zu(:,jlev)+zv(:,jlev)*zv(:,jlev))*0.5 & 04545 & *zp(:)/ga*dsigma(jlev)+zekin(:) 04546 enddo 04547 ! 04548 return 04549 end 04550 subroutine mkekin2(zszp,zsdp,zspp,zekin) 04551 use pumamod 04552 ! 04553 real zszp(NSPP,NLEV),zsdp(NSPP,NLEV),zspp(NSPP) 04554 real zp(NHOR),zekin(NHOR) 04555 ! 04556 real zsd(NESP,NLEV),zsz(NESP,NLEV),zsp(NESP) 04557 real zu(NHOR,NLEV),zv(NHOR,NLEV) 04558 ! 04559 ! some constants 04560 ! 04561 zdelt=delt2/ww ! timestep in s 04562 zcp=gascon/akap ! heat capacity 04563 ! 04564 call mpgallsp(zsd,zsdp,NLEV) 04565 call mpgallsp(zsz,zszp,NLEV) 04566 call mpgallsp(zsp,zspp,NLEV) 04567 do jlev = 1 , NLEV 04568 call dv2uv(zsd(1,jlev),zsz(1,jlev),zu(1,jlev),zv(1,jlev)) 04569 enddo 04570 call sp2fc(zsp,zp) 04571 call fc2gp(zu,NLON,NLPP*NLEV) 04572 call fc2gp(zv,NLON,NLPP*NLEV) 04573 call fc2gp(zp,NLON,NLPP) 04574 ! 04575 zp(:)=psurf*exp(zp(:)) 04576 zekin(:)=0. 04577 do jlev = 1 , NLEV 04578 zu(:,jlev)=cv*zu(:,jlev)*SQRT(rcsq(:)) 04579 zv(:,jlev)=cv*zv(:,jlev)*SQRT(rcsq(:)) 04580 zekin(:)=(zu(:,jlev)*zu(:,jlev)+zv(:,jlev)*zv(:,jlev))*0.5 & 04581 & *zp(:)/ga*dsigma(jlev)+zekin(:) 04582 enddo 04583 ! 04584 return 04585 end 04586 04587 04588 ! ================= 04589 ! SUBROUTINE MKEPOT 04590 ! ================= 04591 04592 subroutine mkepot(zstp,zp,zepot) 04593 use pumamod 04594 ! 04595 real zstp(NSPP,NLEV) 04596 real zp(NHOR),zepot(NHOR) 04597 ! 04598 real zst(NESP,NLEV) 04599 real zt(NHOR,NLEV) 04600 ! 04601 ! some constants 04602 ! 04603 zdelt=delt2/ww ! timestep in s 04604 zcp=gascon/akap ! heat capacity 04605 ! 04606 call mpgallsp(zst,zstp,NLEV) 04607 do jlev = 1 , NLEV 04608 call sp2fc(zst(1,jlev),zt(1,jlev)) 04609 enddo 04610 call fc2gp(zt,NLON,NLPP*NLEV) 04611 ! 04612 zepot(:)=0. 04613 do jlev = 1 , NLEV 04614 zt(:,jlev)=ct*(zt(:,jlev)+t0(jlev)) 04615 zepot(:)=zt(:,jlev)*zcp & 04616 & *zp(:)/ga*dsigma(jlev)+zepot(:) 04617 enddo 04618 ! 04619 return 04620 end 04621 subroutine mkepot2(zstp,zspp,zepot) 04622 use pumamod 04623 ! 04624 real zstp(NSPP,NLEV),zspp(NSPP) 04625 real zp(NHOR),zepot(NHOR) 04626 ! 04627 real zst(NESP,NLEV),zsp(NESP) 04628 real zt(NHOR,NLEV) 04629 ! 04630 ! some constants 04631 ! 04632 zdelt=delt2/ww ! timestep in s 04633 zcp=gascon/akap ! heat capacity 04634 ! 04635 call mpgallsp(zst,zstp,NLEV) 04636 call mpgallsp(zsp,zspp,1) 04637 do jlev = 1 , NLEV 04638 call sp2fc(zst(1,jlev),zt(1,jlev)) 04639 enddo 04640 call sp2fc(zsp,zp) 04641 call fc2gp(zt,NLON,NLPP*NLEV) 04642 call fc2gp(zp,NLON,NLPP) 04643 ! 04644 zp(:)=psurf*exp(zp(:)) 04645 zepot(:)=0. 04646 do jlev = 1 , NLEV 04647 zt(:,jlev)=ct*(zt(:,jlev)+t0(jlev)) 04648 zepot(:)=zt(:,jlev)*zcp & 04649 & *zp(:)/ga*dsigma(jlev)+zepot(:) 04650 enddo 04651 ! 04652 return 04653 end 04654
