 |
PUMA
219
Portable University Model of the Atmosphere
|
|
PUMA
219
Portable University Model of the Atmosphere
|
00001 module pumamod 00002 00003 ! **************************************************************** 00004 ! * Puma Pre Processor * 00005 ! **************************************************************** 00006 ! * E. Kirk & T. Kunz & F. Lunkeit * 00007 ! * Meteorologisches Institut * 00008 ! * Universitaet Hamburg * 00009 ! * Bundesstrasse 55 * 00010 ! * 20146 HAMBURG * 00011 ! * 20-Apr-2009 GERMANY * 00012 ! **************************************************************** 00013 00014 ! **************************************************************** 00015 ! * Insert your own code for modification of initial data to: * 00016 ! * subroutine modify_orography * 00017 ! * subroutine modify_ground_temperature * 00018 ! **************************************************************** 00019 00020 ! **************************************************************** 00021 ! * PUMA in its default setup can run without initial data * 00022 ! * The default setup is an aqua planet with no orography and * 00023 ! * zonally symmetric forcing (newtonian cooling with Tr) * 00024 ! * The atmosphere starts at rest with no horizontal gradients * 00025 ! * * 00026 ! * This preprocessor program performs following tasks: * 00027 ! * 1) Prepare a realistic orography (T21 or T42) * 00028 ! * 2) Enable user modification of this orography * 00029 ! * 3) Enable user modification of the ground temperature field * 00030 ! * 4) Adjust vertical profiles of Restoration Temperature * 00031 ! * 5) Adjust the mean value of surface pressure * 00032 ! * 6) Build an initial Ps field adjusted to orography * 00033 ! * 7) Setup Yoden-profiles * 00034 ! * * 00035 ! * Inputfile: <Naaa_surf_0129.sra> with aaa=032,048,064, ... * 00036 ! * <Naaa_surf_0139.sra> : Ts anomalies * 00037 ! * Outputfiles: <Naaa_surf_0129.sra> : topography [m2/s2] * 00038 ! * <Naaa_surf_0134.sra> : Surface pressure [hPa] * 00039 ! * <Naaa_surf_0121.sra> : Constant part of Tr * 00040 ! * <Naaa_surf_0122.sra> : Variable part of Tr * 00041 ! * <Naaa_surf_0123.sra> : Damping time scales * 00042 ! * * 00043 ! * The outputfiles contain topography, surface pressure, * 00044 ! * the part of Tr, that is constant in time and the part pf Tr * 00045 ! * that can be modulated by an annual cycle. * 00046 ! * * 00047 ! * All files are written formatted, such avoiding the problems * 00048 ! * assigned to big endian and little endian machines * 00049 ! **************************************************************** 00050 00051 00052 ! **************************************************************** 00053 ! * The horizontal resolution of PUMA by the number of latitudes * 00054 ! * nlat is read from file "resolution_namelist" * 00055 ! **************************************************************** 00056 integer :: nlat = 32 00057 00058 ! example values: 32, 48, 64, 128, 192, 256, 512, 1024 00059 ! truncation: T21, T31, T42, T85, T127, T170, T341, T682 00060 00061 00062 ! ***************************************************************** 00063 ! * The number of sigma levels of PUMA are modified after reading * 00064 ! * file <resolution_namelist>. * 00065 ! ***************************************************************** 00066 integer :: nlev = 10 ! Levels 00067 00068 00069 ! *****************************************************************! 00070 ! * Grid related paramters, which are reset after reading the file ! 00071 ! * <resolution_namelist>. All parameters are initialized for the ! 00072 ! * T21 truncation ! 00073 ! *****************************************************************! 00074 integer :: nlon = 64 ! Longitudes = 2 * latitudes 00075 integer :: ntru = 21 ! (nlon-1) / 3 00076 integer :: nlpp = 32 ! Latitudes per process 00077 integer :: nhor = 2048 ! Horizontal part 00078 integer :: nlem = 9 ! Levels - 1 00079 integer :: nlep = 11 ! Levels + 1 00080 integer :: nlsq = 100 ! Levels squared 00081 integer :: ntp1 = 22 ! ntru + 1 00082 integer :: nrsp = 506 ! (ntru+1) * (ntru+2) 00083 integer :: ncsp = 253 ! nrsp / 2 00084 integer :: nspp = 506 ! nodes per process 00085 integer :: nesp = 506 ! number of extended modes 00086 integer :: nzom = 44 ! Number of zonal modes 00087 00088 integer :: nlah = 16 ! Half of latitudes 00089 integer :: nhpp = 16 ! Half latitudes per process 00090 00091 ! **************************************************************** 00092 ! * Don't touch the following parameter definitions ! * 00093 ! **************************************************************** 00094 00095 parameter(AKAP = 0.286) ! Kappa 00096 parameter(ALR = 0.0065) ! Lapse rate 00097 parameter(EZ = 1.632993161855452D0) ! ez = 1 / sqrt(3/8) 00098 parameter(GA = 9.81) ! Gravity 00099 parameter(GASCON = 287.0) ! Gas constant 00100 parameter(PI = 3.141592653589793D0) ! Pi 00101 parameter(TWOPI = PI + PI) ! 2 Pi 00102 parameter(PLARAD_EARTH = 6371000.0) ! Planet radius 00103 parameter(SID_DAY_EARTH= 86164.) ! Siderial day Earth 23h 56m 04s 00104 parameter(PNU = 0.02) ! Time filter 00105 parameter(PNU21 = 1.0 - 2.0*PNU) ! Time filter 2 00106 parameter(PSURF = 101100.0) ! Surface pressure [Pa] 00107 parameter(WW = 0.00007292) ! Rotation speed [1/sec] 00108 parameter(CV = PLARAD_EARTH * WW) ! cv 00109 parameter(CT = CV*CV/GASCON) ! ct 00110 parameter(OSCAR = CV*CV/GA) ! Scale Orography 00111 00112 00113 ! ************* 00114 ! * filenames * 00115 ! ************* 00116 character (256) :: resolution_namelist = "resolution_namelist" 00117 character (256) :: puma_namelist = "puma_namelist" 00118 character (256) :: ppp_puma_txt = "ppp-puma.txt" 00119 00120 ! ***************************************************************** 00121 ! * For multiruns the instance number is appended to the filename * 00122 ! * e.g.: puma_namelist_1 puma_diag_1 etc. for instance # 1 * 00123 ! ***************************************************************** 00124 00125 00126 00127 00128 00129 00130 00131 00132 00133 00134 ! ************************** 00135 ! * Global Integer Scalars * 00136 ! ************************** 00137 00138 integer :: kick = 1 ! kick > 1 initializes eddy generation 00139 integer :: mpstep = 0 ! PUMA 00140 integer :: nafter = 0 ! PUMA 00141 integer :: nwpd = 1 ! PUMA 00142 integer :: ncoeff = 0 ! number of modes to print 00143 integer :: ndel = 6 ! ndel 00144 integer :: ndiag = 12 ! write diagnostics interval 00145 integer :: nkits = 3 ! number of initial timesteps 00146 integer :: nlevt = 9 ! tropospheric levels (set_vertical_grid) 00147 00148 integer :: ngui = 0 ! PUMA variable 00149 integer :: nguidbg = 0 ! PUMA variable 00150 integer :: nqspec = 1 ! PLASIM variable 00151 integer :: nrun = 0 ! PUMA variable 00152 integer :: nstep = 0 ! current timestep 00153 integer :: nstop = 0 ! finishing timestep 00154 integer :: nsync = 0 ! PUMA variable 00155 integer :: ntspd = 24 ! number of timesteps per day 00156 integer :: ncu = 0 ! check unit (debug output) 00157 integer :: noro = 0 ! 1: read orography 00158 ! 2: orography is computed (sine wave) 00159 ! 3: orography is computed (gauss) 00160 00161 integer :: norox = 0 ! x-wavenumber of idealized orography 00162 integer :: nreverse = 0 ! t-gradient reversal in stratosphere; 0:no 1:yes 00163 integer :: ncorrect = 0 ! correct tr due to orography (0:no, 1:yes) 00164 integer :: nsym = 0 ! produces total symmetric initial conditions 00165 00166 integer :: nsrv = 1 ! 1: write gridpoint fields for diagnostics 00167 integer :: lon1oro = 0 ! Define rectangle for anomaly 00168 integer :: lon2oro = 0 ! 00169 integer :: lat1oro = 0 ! 00170 integer :: lat2oro = 0 ! 00171 integer :: ntgr = 0 ! 1: read ground temperature 00172 integer :: lon1tgr = 0 ! Define rectangle for anomaly 00173 integer :: lon2tgr = 0 ! 00174 integer :: lat1tgr = 0 ! 00175 integer :: lat2tgr = 0 ! 00176 integer :: nstrato = 0 ! 1: Torben's stratosphere forcing 00177 integer :: nsponge = 0 ! Switch for sponge layer 00178 integer :: nhelsua = 0 ! 1: Set up Held & Suarez T_R field 00179 ! instead of original PUMA T_R field 00180 ! 2: Set up Held & Suarez T_R field 00181 ! instead of original PUMA T_R field 00182 ! AND use latitudinally varying 00183 ! heating timescale in PUMA (H&Z(94)), 00184 ! irrelevant for PumaPreProcessor (ppp) 00185 ! 3: Use latitudinally varying 00186 ! heating timescale in PUMA (H&Z(94)), 00187 ! irrelevant for PumaPreProcessor (ppp) 00188 integer :: ntestgp = 0 ! switch 0/1 produce a second set of restoration 00189 ! restoration temperatures and damping time scales 00190 integer :: nvg = 0 ! type of vertical grid 00191 ! 0 = linear 00192 ! 1 = Scinocca & Haynes 00193 ! 2 = Polvani & Kushner 00194 00195 integer :: nyoden = 0 ! > 0 Read yoden profile t0(:) and dt(:) 00196 integer :: npackgp = 1 ! used in PUMA 00197 integer :: npacksp = 1 ! used in PUMA 00198 integer :: noutput = 0 ! used in PUMA 00199 integer :: noutsrv = 0 ! used in PUMA 00200 00201 00202 ! These three predifined Yoden profiles may be selected by setting 00203 ! NYODEN= 1, 3 or 5 and nlev=20 00204 00205 real :: t0yod1(20) = (/224.14,213.91,211.36,212.95,217.56 00206 ,224.07,231.24,237.73,243.54,248.87 00207 ,253.41,257.83,261.72,265.60,268.96 00208 ,272.33,275.36,278.34,281.06,283.74/) 00209 real :: dtyod1(20) = (/ 0.00, 0.00, 2.00, 9.19, 19.13 00210 , 28.26, 34.91, 40.37, 45.10, 49.19 00211 , 52.05, 54.69, 56.20, 57.70, 58.31 00212 , 58.91, 59.26, 59.56, 59.76, 59.94/) 00213 real :: t0yod3(20) = (/265.14,254.91,246.36,240.95,237.56 00214 ,234.65,235.24,237.73,243.54,248.87 00215 ,253.41,257.83,261.72,265.60,268.96 00216 ,272.33,275.36,278.34,281.06,283.74/) 00217 real :: dtyod3(20) = (/ 0.00, -3.05,-20.32,-26.19,-28.13 00218 , 28.26, 34.91, 40.37, 45.10, 49.19 00219 , 52.05, 54.69, 56.20, 57.70, 58.31 00220 , 58.91, 59.26, 59.56, 59.76, 59.94/) 00221 real :: t0yod5(20) = (/305.00,298.00,292.00,286.00,280.00 00222 ,274.00,268.00,261.00,256.54,254.87 00223 ,253.41,257.83,261.72,265.60,268.96 00224 ,272.33,275.36,278.34,281.06,283.74/) 00225 real :: dtyod5(20) = (/ 0.00, -3.05,-21.32,-27.19,-29.13 00226 ,-32.26,-34.91,-40.37,-45.10,-49.19 00227 , 52.05, 54.69, 56.20, 57.70, 58.31 00228 , 58.91, 59.26, 59.56, 59.76, 59.94/) 00229 real :: t0yod7(20) = (/265.14,254.91,246.36,240.95,237.56 00230 ,234.65,235.24,237.73,243.54,248.87 00231 ,253.41,257.83,261.72,265.60,268.96 00232 ,272.33,275.36,278.34,281.06,283.74/) 00233 real :: dtyod7(20) = (/ 0.00,-25.05,-64.32,-58.19,-22.13 00234 , 28.26, 34.91, 40.37, 45.10, 49.19 00235 , 52.05, 54.69, 56.20, 57.70, 58.31 00236 , 58.91, 59.26, 59.56, 59.76, 60.00/) 00237 real :: t0yod8(20) = (/265.14,254.91,246.36,240.95,237.56 00238 ,234.65,235.24,237.73,243.54,248.87 00239 ,253.41,257.83,261.72,265.60,268.96 00240 ,272.33,275.36,278.34,281.06,283.74/) 00241 real :: dtyod8(20) = (/ 0.00,-25.05,-64.32,-58.19,-22.13 00242 , 18.26, 24.91, 30.37, 35.10, 39.19 00243 , 42.05, 44.69, 46.20, 47.70, 48.31 00244 , 48.91, 49.26, 49.56, 49.76, 50.00/) 00245 real :: t0yod9(20) = (/265.14,254.91,246.36,240.95,237.56 00246 ,234.65,235.24,237.73,243.54,248.87 00247 ,253.41,257.83,261.72,265.60,268.96 00248 ,272.33,275.36,278.34,281.06,283.74/) 00249 real :: dtyod9(20) = (/ 0.00, 0.00, 0.00, 0.00, 0.00 00250 , 28.26, 34.91, 40.37, 45.10, 49.19 00251 , 52.05, 54.69, 56.20, 57.70, 58.31 00252 , 58.91, 59.26, 59.56, 59.76, 60.00/) 00253 00254 integer :: nfls = 1 ! =1: lower stratospheric forcing is applied, 00255 ! see parameters flsp0, flsdp, flsamp and 00256 ! flsoff below 00257 00258 ! *********************** 00259 ! * Global Real Scalars * 00260 ! *********************** 00261 00262 real :: delt ! 2 pi / ntspd timestep interval 00263 real :: delt2 ! 2 * delt 00264 ! 00265 ! namelist parameter for yoden setup 00266 ! (for arrays to and dt see namelist arrays below) 00267 00268 real :: dtep = 60.0 ! delta T equator <-> pole [K] 00269 real :: dtns = 0.0 ! delta T north <-> south [K] 00270 real :: dtrop = 12000.0 ! Tropopause height [m] 00271 real :: dttrp = 2.0 ! Tropopause smoothing [K] 00272 real :: dtzz = 10.0 ! delta(Theta)/H additional lapserate in 00273 real :: syncstr= 0.0 ! PUMA variable 00274 ! Held & Suarez T_R field 00275 real :: ttp = 200.0 ! Tropopause temperature in 00276 ! Held & Suarez T_R field 00277 real :: plavor= EZ ! planetary vorticity 00278 real :: rotspd = 1.0 ! rotation speed 1.0 = normal Earth rotation 00279 real :: sigmax= 6.0e-7 ! sigma for top half level 00280 real :: tdiss = 0.25 ! diffusion time scale [days] 00281 real :: tac = 0. ! length of annual cycle [days] (0 = no cycle) 00282 real :: pac = 0. ! phase of the annual cycle [days] 00283 real :: tgr = 288.0 ! Ground Temperature in mean profile [K] 00284 real :: oroano= 0.0 ! Orography anomaly in [gpm] 00285 real :: tgrano= 0.0 ! Ground temperature anomaly in [K] 00286 !! 00287 real :: alrs = -0.0000! stratospheric lapse rate [K/m] 00288 real :: horo = 0.0 ! height of the idealized orography [m] 00289 real :: orofac= 1.0 ! factor to scale the orograpy 00290 real :: ttropo= 250.0 ! temp. at tropopause 00291 !! 00292 real :: dorox = 0.0 ! lon. base point for gauss mountain [dec] 00293 real :: doroy = 0.0 ! lat. base point for gauss mountain [dec] 00294 real :: doroxs = 0.0 ! gauss mountain scale in lon.-dir. [dec] 00295 real :: doroys = 0.0 ! gauss mountain scale in lat.-dir. [dec] 00296 00297 real :: tauta = 40.0 ! heating timescale far from surface [days] 00298 real :: tauts = 4.0 ! heating timescale close to surface [days] 00299 real :: sid_day= SID_DAY_EARTH ! siderial day [sec] 00300 00301 00302 00303 real :: zsigb = 0.7 ! sigma_b for Held&Suarez frict. and heat. 00304 ! time scale 00305 00306 00307 00308 ! * parameter for stratospheric forcing * 00309 00310 real :: alrpv = 0.002 ! Vertical lapse rate of stratospheric 00311 ! polar vortex forcing restoration 00312 ! temperature field. It corresponds to 00313 ! 'gamma' in Polvani & Kushner (2002). 00314 ! But alrpv is in [K/m], thus, 00315 ! alrpv=0.002 corresponds to 'gamma=2' 00316 ! in P&K (2002). 00317 real :: radpv = 50. ! Radius of stratospheric polar vortex 00318 ! forcing [deg latitude] 00319 real :: edgepv = 10. ! Width of edge of stratospheric polar 00320 ! vortex forcing [deg latitude]. 00321 ! If edgepv=0., then no polar vortex is 00322 ! set up. 00323 real :: pmaxpv = 10000. ! Lower boundary of stratospheric polar 00324 ! vortex forcing, max. pressure [Pa]. 00325 ! If pmaxpv=0. then pmaxpv is set to the 00326 ! pressure at tropopause height, specified 00327 ! by dtrop, according to the US standard 00328 ! atmosphere (USSA) vertical profile, used 00329 ! for contruction of the restoration 00330 ! temperature field. With the standard 00331 ! setting (dtrop=11000m, ALR=0.0065K/m, 00332 ! tgr=288.15K) this gives pmaxpv=22632Pa. 00333 ! ++++NOTE: pmaxpv should be within the 00334 ! second USSA layer, that is the interval 00335 ! from 54749Pa to 22632Pa (in case of 00336 ! standard parameter setting).++++ 00337 00338 real :: flsp0 = 10000. ! pressure of max. lower stratospheric forcing 00339 ! field [Pa] 00340 real :: flsdp = 9000. ! half of vertical extension of lower 00341 ! stratospheric forcing field [Pa] 00342 real :: flsamp = 15. ! amplitude of lower stratospheric forcing 00343 ! field [K] 00344 real :: flsoff = -5. ! offset of lower stratospheric forcing 00345 ! field [K] 00346 00347 00348 ! ******************* 00349 ! * Namelist Arrays * 00350 ! ******************* 00351 ! workaround for older fortran versions, where allocatable arrays 00352 ! are not allowed in namelists 00353 00354 00355 integer, parameter :: MAXLEV = 100 00356 integer, parameter :: MAXSELZW = 42 00357 integer, parameter :: MAXSELSP = ((MAXSELZW+1) * (MAXSELZW+2)) / 2 00358 integer :: nselect(0:MAXSELZW) = 1 ! NSELECT can be used up tp T42 00359 integer :: nspecsel(MAXSELSP) = 1 00360 integer :: ndl(MAXLEV) = 0 00361 real :: restim(MAXLEV) = 15.0 00362 real :: sigmah(MAXLEV) = 0.0 00363 real :: t0k(MAXLEV) = 250.0 00364 real :: t0(MAXLEV) = 250.0 00365 real :: tfrc(MAXLEV) = 0.0 00366 real :: dt(MAXLEV) = 0.0 00367 00368 ! ************************** 00369 ! * Global Spectral Arrays * 00370 ! ************************** 00371 00372 real, allocatable :: sor(:) ! Spectral Orography 00373 real, allocatable :: ssp(:) ! Spectral surface pressure 00374 real, allocatable :: stg(:) ! Spectral ground temperature 00375 real, allocatable :: sep(:) ! Spectral equator-pole gradient 00376 real, allocatable :: sns(:) ! Spectral north-south gradient 00377 real, allocatable :: spnorm(:) ! ECHAM -> PUMA normalization factors 00378 real, allocatable :: sr1(:,:) ! Constant part of Tr 00379 real, allocatable :: sr2(:,:) ! Variable part of Tr 00380 00381 ! *************************** 00382 ! * Global Gridpoint Arrays * 00383 ! *************************** 00384 00385 real, allocatable :: gor(:,:) ! Orography 00386 real, allocatable :: gsp(:,:) ! Surface pressure 00387 real, allocatable :: gtg(:,:) ! Ground temperature 00388 real, allocatable :: gan(:,:) ! Ground temperature anomaly 00389 real, allocatable :: gep(:,:) ! Equator-pole gradient 00390 real, allocatable :: gns(:,:) ! North-south gradient 00391 real, allocatable :: gtc(:,:,:) ! Restoration Temperature 00392 real, allocatable :: gtv(:,:,:) ! Restoration Temperature NS mode 00393 real, allocatable :: gtdamp(:,:,:)! Reciprocal of damping time scale for T 00394 real, allocatable :: gaf(:,:,:) ! Anomaly factors 00395 real, allocatable :: gra(:) ! Gradient factors 00396 00397 ! ******************* 00398 ! * Latitude Arrays * 00399 ! ******************* 00400 00401 character (3), allocatable :: chlat(:) ! label for latitudes 00402 00403 real (kind=8), allocatable :: sid(:) ! sin(phi) 00404 real (kind=8), allocatable :: gwd(:) ! Gaussian weight 00405 real (kind=8), allocatable :: csq(:) ! cos(phi)^2 00406 real (kind=8), allocatable :: dla(:) ! phi 00407 00408 00409 ! **************** 00410 ! * Level Arrays * 00411 ! **************** 00412 00413 real, allocatable :: t0d(:) ! vertical t0k gradient 00414 real, allocatable :: dsigma(:) 00415 real, allocatable :: rdsig(:) 00416 real, allocatable :: sigma(:) 00417 real, allocatable :: tkp(:) 00418 00419 ! ********************** 00420 ! * Dummy declarations * 00421 ! ********************** 00422 00423 real :: gp(2) 00424 real :: sp(2) 00425 real :: gpj(2) 00426 real :: gu(2,2) 00427 real :: gv(2,2) 00428 real :: sd(2,2) 00429 real :: st(2,2) 00430 real :: sq(2,2) 00431 real :: sz(2,2) 00432 real :: gd(2,2) 00433 real :: gq(2,2) 00434 real :: gt(2,2) 00435 real :: gz(2,2) 00436 00437 end module pumamod 00438 00439 program ppp 00440 00441 ! open file ppp-puma interface information 00442 if (mypid == NROOT) then 00443 open(95,file="ppp-puma.txt",form='formatted') 00444 endif 00445 00446 call resolution 00447 call prolog 00448 call gridpoint 00449 00450 ! close file ppp-puma interface information 00451 if (mypid == NROOT) then 00452 close(95) 00453 endif 00454 00455 stop 00456 end 00457 00458 ! ****************************** 00459 ! * SUBROUTINE ALLOCATE_ARRAYS * 00460 ! ****************************** 00461 00462 subroutine allocate_arrays 00463 use pumamod 00464 00465 !--- Global Spectral Arrays 00466 allocate(sor(nesp)) ; sor(:) = 0.0 ! Spectral Orography 00467 allocate(ssp(nesp)) ; ssp(:) = 0.0 ! Spectral surface pressure 00468 allocate(stg(nesp)) ; stg(:) = 0.0 ! Spectral ground temperature 00469 allocate(sep(nesp)) ; sep(:) = 0.0 ! Spectral equator-pole gradient 00470 allocate(sns(nesp)) ; sns(:) = 0.0 ! Spectral north-south gradient 00471 allocate(spnorm(nesp)) ; spnorm(:) = 0.0 ! ECHAM -> PUMA normalization 00472 allocate(sr1(nesp,nlev)) ; sr1(:,:) = 0.0 ! Constant part of Tr 00473 allocate(sr2(nesp,nlev)) ; sr2(:,:) = 0.0 ! Variable part of Tr 00474 00475 !--- Global Gridpoint Arrays 00476 allocate(gor(nlon,nlat)) ; gor(:,:) = 0.0 ! Orography 00477 allocate(gsp(nlon,nlat)) ; gsp(:,:) = 0.0 ! Surface pressure 00478 allocate(gtg(nlon,nlat)) ; gtg(:,:) = 0.0 ! Ground temperature 00479 allocate(gan(nlon,nlat)) ; gan(:,:) = 0.0 ! Ground temperature anomaly 00480 allocate(gep(nlon,nlat)) ; gep(:,:) = 0.0 ! Equator-pole gradient 00481 allocate(gns(nlon,nlat)) ; gns(:,:) = 0.0 ! North-south gradient 00482 allocate(gtc(nlon,nlat,nlev)) ; gtc(:,:,:) = 0.0 ! Restoration Temperature 00483 allocate(gtv(nlon,nlat,nlev)) ; gtv(:,:,:) = 0.0 ! Rest. Temperature NS mode 00484 allocate(gtdamp(nlon,nlat,nlev)); gtdamp(:,:,:) = 0.0 ! Reciprocal of damping time scale for T 00485 allocate(gaf(nlon,nlat,nlev)) ; gaf(:,:,:) = 0.0 ! Anomaly factors 00486 00487 !--- Latitude Arrays 00488 allocate(chlat(nlat)) ; chlat(:) = " " ! label for latitudes 00489 allocate(sid(nlat)) ; sid(:) = 0.0 ! sin(phi) 00490 allocate(gwd(nlat)) ; gwd(:) = 0.0 ! Gaussian weight 00491 allocate(csq(nlat)) ; csq(:) = 0.0 ! cos(phi)^2 00492 allocate(dla(nlat)) ; dla(:) = 0.0 ! phi 00493 00494 !--- Level Arrays 00495 allocate(gra(nlev)) ; gra(:) = 0.0 ! Gradient factors 00496 allocate(t0d(nlev)) ; t0d(:) = 0.0 ! vertical t0k gradient 00497 allocate(dsigma(nlev)) ; dsigma(:) = 0.0 00498 allocate(rdsig(nlev)) ; rdsig(:) = 0.0 00499 allocate(sigma(nlev)) ; sigma(:) = 0.0 00500 allocate(tkp(nlev)) ; tkp(:) = 0.0 00501 00502 return 00503 end subroutine allocate_arrays 00504 00505 ! =========================== 00506 ! SUBROUTINE MODIFY_OROGRAPHY 00507 ! =========================== 00508 00509 subroutine modify_orography(por) 00510 use pumamod 00511 real :: por(nlon,nlat) 00512 00513 ! Array <por> contains the orography on a Gaussian grid 00514 ! the units are that of a geopotential [m2/s2] (gpm * g) 00515 ! You may modify the orography here with your own code 00516 ! The new orography is spectrally fitted after this routine 00517 ! A gridpoint representation is written in service format 00518 ! to file <puma_oro_tnn.srv> with nn = spectral truncation. 00519 00520 if (noro == 2) call mkoro(por) 00521 if (noro == 3) call mkorog(por) 00522 00523 ! Rectangular anomaly 00524 00525 ! if (oroano /= 0.0 .and. lat1oro > 0 .and. lat2oro <= nlat & 00526 ! .and. lon1oro > 0 .and. lon2oro <= nlon & 00527 ! .and. lat1oro < lat2oro .and. lon1oro < lon2oro) then 00528 ! do jlat = lat1oro , lat2oro 00529 ! do jlon = lon1oro , lon2oro 00530 ! por(jlon,jlat) = por(jlon,jlat) + oroano * GA 00531 ! if (por(jlon,jlat) < 0.0) por(jlon,jlat) = 0.0 00532 ! enddo 00533 ! enddo 00534 ! endif 00535 00536 ! Elliptic anomaly 00537 00538 if (oroano /= 0.0 .and. lat1oro > 0 .and. lat2oro <= nlat & 00539 .and. lon1oro > 0 .and. lon2oro <= nlon & 00540 .and. lat1oro < lat2oro .and. lon1oro < lon2oro) then 00541 x0 = (lon1oro + lon2oro) * 0.5 00542 y0 = (lat1oro + lat2oro) * 0.5 00543 xf = PI / (lon2oro - lon1oro) 00544 yf = PI / (lat2oro - lat1oro) 00545 do jlat = lat1oro , lat2oro 00546 yb = (jlat - y0) * yf 00547 do jlon = lon1oro , lon2oro 00548 xa = (jlon - x0) * xf 00549 cx = cos(xa) 00550 cy = cos(yb) 00551 if (cx > 0.0 .and. cy > 0.0) then 00552 por(jlon,jlat) = por(jlon,jlat) + oroano * GA * cx * cy 00553 endif 00554 if (por(jlon,jlat) < 0.0) por(jlon,jlat) = 0.0 00555 enddo 00556 enddo 00557 endif 00558 00559 return 00560 end 00561 00562 ! ==================================== 00563 ! SUBROUTINE MODIFY_GROUND_TEMPERATURE 00564 ! ==================================== 00565 00566 subroutine modify_ground_temperature(ptgr) 00567 use pumamod 00568 real :: ptgr(nlon,nlat) 00569 00570 ! Array <ptgr> contains the ground temperature on a Gaussian grid. 00571 ! The units are [K] 00572 ! You may modify the ground temperature here with your own code. 00573 ! The new ground temperature is used to construct the temperature 00574 ! profile of the restoration temperature on each grid point. 00575 ! A gridpoint representation is written in service format 00576 ! to file <puma_gtgr_tnn.srv> with nn = spectral truncation. 00577 00578 if (tgrano /= 0.0 .and. lat1tgr > 0 .and. lat2tgr <= nlat & 00579 .and. lon1tgr > 0 .and. lon2tgr <= nlon & 00580 .and. lat1tgr < lat2tgr .and. lon1tgr < lon2tgr) then 00581 do jlat = lat1tgr , lat2tgr 00582 do jlon = lon1tgr , lon2tgr 00583 gan(jlon,jlat) = gan(jlon,jlat) + tgrano 00584 ptgr(jlon,jlat) = ptgr(jlon,jlat) + tgrano 00585 if (ptgr(jlon,jlat) < 0.0) ptgr(jlon,jlat) = 0.0 00586 enddo 00587 enddo 00588 endif 00589 00590 return 00591 end 00592 00593 00594 ! ================ 00595 ! SUBROUTINE MKORO 00596 ! ================ 00597 00598 subroutine mkoro(por) 00599 use pumamod 00600 ! 00601 real por(nlon,nlat) 00602 ! 00603 zscale=horo*GA*0.5 00604 ! 00605 por(:,:)=0. 00606 ! 00607 do jlat=1,nlat/2 00608 zlat=dla(jlat) 00609 zfacy=(sin(2.*zlat))**2 00610 do jlon=1,nlon 00611 if(norox > 0) then 00612 zfacx=(1.+cos(norox*real(jlon)*TWOPI/real(nlon))) 00613 else 00614 zfacx=1. 00615 endif 00616 por(jlon,jlat)=zscale*zfacx*zfacy 00617 enddo 00618 enddo 00619 ! 00620 if(nsym == 1) then 00621 do jlat=1,nlat/2 00622 j2=nlat+1-jlat 00623 por(:,j2)=por(:,jlat) 00624 enddo 00625 endif 00626 ! 00627 return 00628 end 00629 00630 ! ================= 00631 ! SUBROUTINE MKOROG 00632 ! ================= 00633 00634 subroutine mkorog(por) 00635 use pumamod 00636 ! 00637 real por(nlon,nlat) 00638 ! 00639 zscale=horo*GA 00640 ! 00641 zlon0=dorox*PI/180. 00642 zlat0=doroy*PI/180. 00643 zlons=(180./(doroxs*PI))**2 00644 zlats=(180./(doroys*PI))**2 00645 ! 00646 do jlat=1,nlat 00647 zlat=dla(jlat) 00648 !! zcos2=cos(zlat)**2 00649 zcos2=1. 00650 zdlat2=(zlat-zlat0)**2 00651 do jlon=1,nlon 00652 zlon=TWOPI*real(jlon-1)/real(nlon) 00653 zdlon=abs(zlon-zlon0) 00654 if(zdlon > PI) zdlon=TWOPI-zdlon 00655 zdlon2=zdlon**2 00656 por(jlon,jlat)=zscale*EXP(-zlons*zcos2*zdlon2-zlats*zdlat2) 00657 enddo 00658 enddo 00659 ! 00660 if(nsym == 1 .and. zlat0 > 0.) then 00661 do jlat=1,nlat/2 00662 j2=nlat+1-jlat 00663 por(:,j2)=por(:,jlat) 00664 enddo 00665 elseif(nsym == 1 .and. zlat0 < 0.) then 00666 do jlat=1,nlat/2 00667 j2=nlat+1-jlat 00668 por(:,jlat)=por(:,j2) 00669 enddo 00670 endif 00671 ! 00672 return 00673 end 00674 00675 ! ================= 00676 ! SUBROUTINE PROLOG 00677 ! ================= 00678 00679 subroutine prolog 00680 use pumamod 00681 00682 00683 00684 call allocate_arrays 00685 00686 call printparameter 00687 call inigau(nlat,sid,gwd) 00688 call inilat 00689 call legpri 00690 call readnl 00691 call initpm 00692 call legini(nlat,nlpp,nesp,nlev,plavor,sid,gwd) 00693 call initfd 00694 00695 00696 if (mypid == NROOT) then 00697 call ppp_write_i('NLAT',1,nlat) 00698 call ppp_write_i('NLEV',1,nlev) 00699 call ppp_write_i('NHELSUA',1,nhelsua) 00700 call ppp_write_r('SIGMH',NLEV,sigmah) 00701 endif 00702 00703 return 00704 end 00705 00706 ! ================= 00707 ! SUBROUTINE INITFD 00708 ! ================= 00709 00710 subroutine initfd 00711 use pumamod 00712 00713 dimension zrmean(nlev) 00714 00715 zfmode0 = sqrt(2.0) 00716 zfmode1 = 1.0 / sqrt(6.0) 00717 zfmode2 = -2.0 / 3.0 * sqrt(0.4) 00718 00719 stg(1) = zfmode0 * tgr ! Ground temperature [K] 00720 sns(3) = zfmode1 * dtns ! North-South gradient [K] 00721 sep(5) = zfmode2 * dtep ! Equator-Pole gradient [K] 00722 00723 ! Find sigma at dtrop 00724 00725 zttrop = tgr - dtrop * ALR 00726 ztps = (zttrop/tgr)**(GA/(ALR*GASCON)) 00727 00728 ! The North-South and Equator-Pole gradients are defined on z=0. 00729 ! gra() modifies the gradient from full at z=0 to zero at tropopause 00730 ! PUMA aquaplanet compatibility mode (sine function used) 00731 00732 do jlev = 1 , nlev 00733 gra(jlev) = sin(0.5 * PI * (sigma(jlev) - ztps) / (1.0-ztps)) 00734 if (gra(jlev) < 0.0) gra(jlev) = 0.0 00735 enddo 00736 00737 return 00738 end 00739 00740 ! ========================== 00741 ! SUBROUTINE ANOMALY_FACTORS 00742 ! ========================== 00743 00744 subroutine anomaly_factors 00745 use pumamod 00746 00747 do jlat = 1 , nlat 00748 do jlon = 1 , nlon 00749 00750 ! Find sigma at dtrop 00751 00752 zttrop = tgr - dtrop * ALR 00753 ztps = (zttrop/gtg(jlon,jlat))**(GA/(ALR*GASCON)) 00754 00755 ! The North-South and Equator-Pole gradients are defined on z=0. 00756 ! gaf() modifies the gradient from full at z=0 to zero at tropopause 00757 ! PUMA aquaplanet compatibility mode (sine function used) 00758 00759 if(nreverse == 0) then 00760 do jlev = 1 , nlev 00761 gaf(jlon,jlat,jlev) = sin(0.5 * PI * (sigma(jlev) - ztps) / (1.0-ztps)) 00762 if (gaf(jlon,jlat,jlev) < 0.0) gaf(jlon,jlat,jlev) = 0.0 00763 enddo 00764 else 00765 do jlev = 1 , nlev 00766 gaf(jlon,jlat,jlev) = sin(0.5 * PI * (sigma(jlev) - ztps) / (1.0-ztps)) 00767 if (sigma(jlev) < ztps) & 00768 & gaf(jlon,jlat,jlev) = sin(PI*(sigma(jlev)-ztps)/(ztps-sigma(1))) 00769 enddo 00770 endif 00771 00772 enddo 00773 enddo 00774 00775 return 00776 end 00777 00778 00779 ! ======================== 00780 ! SUBROUTINE READ_GAN_GRID 00781 ! ======================== 00782 00783 subroutine read_gan_grid(kread) 00784 use pumamod 00785 00786 logical :: lexist 00787 integer :: ihead(8) 00788 character(20) :: ynumber 00789 character(256) :: yfilename 00790 00791 ! Read temperature anomalies for PUMA 00792 ! Use formatted Service style 00793 00794 kread = 0 00795 00796 if (nlat < 1000) then 00797 write(ynumber,'(I3.3)') nlat 00798 else 00799 write(ynumber,'(I4.4)') nlat 00800 endif 00801 00802 yfilename = "N" // trim(adjustl(ynumber)) // "_surf_0139.sra" 00803 inquire(file=yfilename,exist=lexist) 00804 if (lexist) then 00805 write(*,*) ' Reading anomaly temperature from file <',trim(yfilename),'>' 00806 open (65,file=yfilename,form='formatted') 00807 read (65,*) ihead(:) 00808 read (65,*) gan(:,:) 00809 close(65) 00810 kread = 1 00811 else 00812 write(*,*) ' No anomaly temperature file' 00813 endif 00814 return 00815 end subroutine read_gan_grid 00816 00817 00818 ! ======================== 00819 ! SUBROUTINE READ_ORO_GRID 00820 ! ======================== 00821 00822 subroutine read_oro_grid(kread) 00823 use pumamod 00824 00825 logical :: lexist 00826 integer :: ihead(8) 00827 character(20) :: ynumber 00828 character(256) :: yfilename 00829 00830 ! Read orography for PUMA 00831 ! Use formatted Service style 00832 00833 kread = 0 00834 00835 if (nlat < 1000) then 00836 write(ynumber,'(I3.3)') nlat 00837 else 00838 write(ynumber,'(I4.4)') nlat 00839 endif 00840 00841 yfilename = "N" // trim(adjustl(ynumber)) // "_surf_0129.sra" 00842 inquire(file=yfilename,exist=lexist) 00843 if (lexist) then 00844 write(*,*) ' Reading orography from file <',trim(yfilename),'>' 00845 open (65,file=yfilename,form='formatted') 00846 read (65,*) ihead(:) 00847 read (65,*) gor(:,:) 00848 close(65) 00849 kread = 1 00850 else 00851 gor(:,:) = 0.0 00852 write(*,*) ' No orography file - starting with zero orography' 00853 endif 00854 return 00855 end subroutine read_oro_grid 00856 00857 00858 ! ==================== 00859 ! SUBROUTINE WRITE_ORO 00860 ! ==================== 00861 00862 subroutine write_oro 00863 use pumamod 00864 00865 dimension itime(8) 00866 dimension ihead(8) 00867 00868 character(20) :: ynumber 00869 character(256) :: yfilename 00870 00871 call date_and_time(values=itime) 00872 00873 ! Write orography for PUMA 00874 ! Use formatted Service style 00875 00876 if (nlat < 1000) then 00877 write(ynumber,'(I3.3)') nlat 00878 else 00879 write(ynumber,'(I4.4)') nlat 00880 endif 00881 00882 yfilename = "N" // trim(adjustl(ynumber)) // "_surf_0129.sra" 00883 open(60,file=yfilename,form='formatted') 00884 ihead(1) = 129 ! code for orography 00885 ihead(2) = 0 ! level 00886 ihead(3) = itime(1) * 10000 + itime(2) * 100 + itime(3) ! YYYYMMDD 00887 ihead(4) = itime(5) * 100 + itime(6) ! HHMM 00888 ihead(5) = nlon ! 1. dimension 00889 ihead(6) = nlat ! 2. dimension 00890 ihead(7) = 1 00891 ihead(8) = 0 00892 00893 write (60,'(8I10)') ihead(:) 00894 write (60,'(8F10.3)') gor(:,:) 00895 00896 close(60) 00897 00898 return 00899 end subroutine write_oro 00900 00901 00902 ! =================== 00903 ! SUBROUTINE WRITE_PS 00904 ! =================== 00905 00906 subroutine write_ps 00907 use pumamod 00908 00909 dimension itime(8) 00910 dimension ihead(8) 00911 00912 character(20) :: ynumber 00913 character(256) :: yfilename 00914 00915 real :: zpres(nlon,nlat) 00916 00917 call date_and_time(values=itime) 00918 00919 ! Write surface pressure for PUMA 00920 ! Use formatted Service style 00921 00922 if (nlat < 1000) then 00923 write(ynumber,'(I3.3)') nlat 00924 else 00925 write(ynumber,'(I4.4)') nlat 00926 endif 00927 yfilename = "N" // trim(adjustl(ynumber)) // "_surf_0134.sra" 00928 open(60,file=yfilename,form='formatted') 00929 ihead(1) = 134 ! code for surface pressure [hPa] 00930 ihead(2) = 0 ! level 00931 ihead(3) = itime(1) * 10000 + itime(2) * 100 + itime(3) ! YYYYMMDD 00932 ihead(4) = itime(5) * 100 + itime(6) ! HHMM 00933 ihead(5) = nlon ! 1. dimension 00934 ihead(6) = nlat ! 2. dimension 00935 ihead(7) = 1 00936 ihead(8) = 0 00937 00938 zpres(:,:) = 0.01 * PSURF * exp(gsp(:,:)) ! Store as [hPa] 00939 00940 write (60,'(8I10)') ihead(:) 00941 write (60,'(8F10.4)') zpres(:,:) 00942 00943 close(60) 00944 00945 return 00946 end subroutine write_ps 00947 00948 00949 ! ==================== 00950 ! SUBROUTINE WRITE_GTC 00951 ! ==================== 00952 00953 subroutine write_gtc 00954 use pumamod 00955 00956 dimension itime(8) 00957 dimension ihead(8) 00958 00959 character(20) :: ynumber 00960 character(256) :: yfilename 00961 00962 call date_and_time(values=itime) 00963 00964 ! Write constant part of Tr 00965 ! Use formatted Service style 00966 00967 if (nlat < 1000) then 00968 write(ynumber,'(I3.3)') nlat 00969 else 00970 write(ynumber,'(I4.4)') nlat 00971 endif 00972 yfilename = "N" // trim(adjustl(ynumber)) // "_surf_0121.sra" 00973 00974 open(60,file=yfilename,form='formatted') 00975 ihead(1) = 121 ! code for Tr const 00976 ihead(2) = 0 ! level 00977 ihead(3) = itime(1) * 10000 + itime(2) * 100 + itime(3) ! YYYYMMDD 00978 ihead(4) = itime(5) * 100 + itime(6) ! HHMM 00979 ihead(5) = nlon ! 1. dimension 00980 ihead(6) = nlat ! 2. dimension 00981 ihead(7) = 1 00982 ihead(8) = 0 00983 00984 do jlev = 1 , nlev 00985 ihead(2) = jlev 00986 write (60,'(8I10)') ihead(:) 00987 write (60,'(8F10.4)') gtc(:,:,jlev) 00988 enddo 00989 00990 close(60) 00991 00992 return 00993 end subroutine write_gtc 00994 00995 00996 00997 ! ==================== 00998 ! SUBROUTINE WRITE_GTV 00999 ! ==================== 01000 01001 subroutine write_gtv 01002 use pumamod 01003 01004 dimension itime(8) 01005 dimension ihead(8) 01006 01007 character(20) :: ynumber 01008 character(256) :: yfilename 01009 01010 call date_and_time(values=itime) 01011 01012 ! Write variable part of Tr 01013 ! Use formatted Service style 01014 01015 if (nlat < 1000) then 01016 write(ynumber,'(I3.3)') nlat 01017 else 01018 write(ynumber,'(I4.4)') nlat 01019 endif 01020 yfilename = "N" // trim(adjustl(ynumber)) // "_surf_0122.sra" 01021 open(60,file=yfilename,form='formatted') 01022 ihead(1) = 122 ! code for Tr variable 01023 ihead(2) = 0 ! level 01024 ihead(3) = itime(1) * 10000 + itime(2) * 100 + itime(3) ! YYYYMMDD 01025 ihead(4) = itime(5) * 100 + itime(6) ! HHMM 01026 ihead(5) = nlon ! 1. dimension 01027 ihead(6) = nlat ! 2. dimension 01028 ihead(7) = 1 01029 ihead(8) = 0 01030 01031 do jlev = 1 , nlev 01032 ihead(2) = jlev 01033 write (60,'(8I10)') ihead(:) 01034 write (60,'(8F10.4)') gtv(:,:,jlev) 01035 enddo 01036 01037 close(60) 01038 01039 return 01040 end subroutine write_gtv 01041 01042 01043 ! ======================== 01044 ! SUBROUTINE WRITE_VARGP2D 01045 ! ======================== 01046 01047 subroutine write_vargp2D(zgp,kcode) 01048 use pumamod 01049 01050 dimension itime(8) 01051 dimension ihead(8) 01052 01053 character(20) :: ynumber 01054 character(256) :: yfilename 01055 real :: zgp(nlon,nlat) 01056 01057 01058 call date_and_time(values=itime) 01059 01060 ! produce file name to be written 01061 if (NLAT < 1000) then 01062 write(yfilename,'("N",I3.3,"_surf_",I4.4,".sra")') NLAT,kcode 01063 else 01064 write(yfilename,'("N",I4.4,"_surf_",I4.4,".sra")') NLAT,kcode 01065 endif 01066 01067 01068 01069 open(60,file=yfilename,form='formatted') 01070 ihead(1) = kcode ! code for reciprocal of damping time scale 01071 ihead(2) = 0 ! level 01072 ihead(3) = itime(1) * 10000 + itime(2) * 100 + itime(3) ! YYYYMMDD 01073 ihead(4) = itime(5) * 100 + itime(6) ! HHMM 01074 ihead(5) = nlon ! 1. dimension 01075 ihead(6) = nlat ! 2. dimension 01076 ihead(7) = 1 01077 ihead(8) = 0 01078 01079 select case(kcode) 01080 case(129,134) 01081 ihead(2) = 0 01082 write (60,'(8I10)') ihead(:) 01083 write (60,'(8(X,E16.10))') zgp(:,:) 01084 end select 01085 01086 close(60) 01087 01088 return 01089 end subroutine write_vargp2D 01090 01091 ! ======================== 01092 ! SUBROUTINE WRITE_VARGP3D 01093 ! ======================== 01094 01095 subroutine write_vargp3D(zgp,kcode,klev) 01096 use pumamod 01097 01098 dimension itime(8) 01099 dimension ihead(8) 01100 01101 character(20) :: ynumber 01102 character(256) :: yfilename 01103 real :: zgp(nlon,nlat,klev) 01104 01105 01106 call date_and_time(values=itime) 01107 01108 ! produce file name to be written 01109 if (NLAT < 1000) then 01110 write(yfilename,'("N",I3.3,"_surf_",I4.4,".sra")') NLAT,kcode 01111 else 01112 write(yfilename,'("N",I4.4,"_surf_",I4.4,".sra")') NLAT,kcode 01113 endif 01114 01115 01116 01117 open(60,file=yfilename,form='formatted') 01118 ihead(1) = kcode ! code for reciprocal of damping time scale 01119 ihead(2) = 0 ! level 01120 ihead(3) = itime(1) * 10000 + itime(2) * 100 + itime(3) ! YYYYMMDD 01121 ihead(4) = itime(5) * 100 + itime(6) ! HHMM 01122 ihead(5) = nlon ! 1. dimension 01123 ihead(6) = nlat ! 2. dimension 01124 ihead(7) = 1 01125 ihead(8) = 0 01126 01127 select case(kcode) 01128 case(121,122,123,124,125,126) 01129 do jlev = 1 , klev 01130 ihead(2) = jlev 01131 write (60,'(8I10)') ihead(:) 01132 write (60,'(8(X,E16.10))') zgp(:,:,jlev) 01133 enddo 01134 end select 01135 01136 close(60) 01137 01138 return 01139 end subroutine write_vargp3D 01140 01141 01142 ! ==================== 01143 ! SUBROUTINE GRIDPOINT 01144 ! ==================== 01145 01146 subroutine gridpoint 01147 use pumamod 01148 !! 01149 dimension ihead(8) 01150 dimension zprof(nlev) 01151 dimension zzprof(nlon,nlat,nlev) 01152 character(256) :: yfilename,yfohead,yfodata,ymessage 01153 logical :: exist 01154 01155 ! Read orography 01156 01157 gor(:,:) = 0.0 01158 if (noro > 0) then 01159 call read_oro_grid(iread) 01160 endif 01161 01162 ! Read ground anomaly temperature 01163 01164 call read_gan_grid(iread) 01165 01166 call sp2fc(sns,gns) 01167 call sp2fc(sep,gep) 01168 call sp2fc(stg,gtg) 01169 01170 call alt2reg(gns,1) 01171 call alt2reg(gep,1) 01172 call alt2reg(gtg,1) 01173 01174 01175 call fc2gp(gns,nlon,nlat) 01176 call fc2gp(gep,nlon,nlat) 01177 call fc2gp(gtg,nlon,nlat) 01178 01179 gtg(:,:) = gtg(:,:) + gan(:,:) 01180 01181 call modify_orography(gor) ! User interface 01182 01183 if (nyoden /= 0) then 01184 call yoden 01185 else 01186 ! compute ground temperature on orography surface 01187 01188 gtg(:,:) = gtg(:,:) - ALR * gor(:,:) / GA ! [gpm] 01189 01190 call modify_ground_temperature(gtg) ! User interface 01191 01192 call anomaly_factors ! Compute factors for NS & EP 01193 01194 ! Compute vertical profile for each column 01195 01196 do jlat = 1 , nlat 01197 do jlon = 1 , nlon 01198 call tprofile(gtg(jlon,jlat),zprof,gor(jlon,jlat)/GA) 01199 gtc(jlon,jlat,:) = zprof(:) 01200 enddo 01201 enddo 01202 01203 ! Modify Restoration Temperature with EP mode 01204 01205 do jlev = 1 , nlev 01206 gtc(:,:,jlev) = gtc(:,:,jlev) + gaf(:,:,jlev) * gep(:,:) 01207 enddo 01208 01209 ! Compute vertical profile for each column including variable NS mode 01210 01211 do jlev = 1 , nlev 01212 gtv(:,:,jlev) = gaf(:,:,jlev) * gns(:,:) 01213 enddo 01214 01215 endif ! (nyoden == 1) 01216 01217 ! Initialize surface pressure (LnPs) 01218 01219 do jlat = 1 , nlat 01220 do jlon = 1 , nlon 01221 gsp(jlon,jlat) = -gor(jlon,jlat) / (GASCON*tgr) 01222 enddo 01223 enddo 01224 01225 ! if (nhelsua == 1 .or. nhelsua == 2) then 01226 if (nhelsua > 0) then 01227 call heldsuarez 01228 gor(:,:) = 0.0 01229 gsp(:,:) = 0.0 01230 endif 01231 01232 01233 if (nstrato == 1) then 01234 call setzt2 ! Torben's forcing initialisation 01235 endif 01236 01237 call write_oro 01238 call write_ps 01239 call write_gtc 01240 call write_gtv 01241 call write_vargp3D(gtdamp,123,nlev) 01242 if (ntestgp == 1) then 01243 call write_vargp3D(gtc,124,nlev) 01244 call write_vargp3D(gtv,125,nlev) 01245 call write_vargp3D(gtdamp,126,nlev) 01246 endif 01247 01248 call printprofile 01249 01250 return 01251 end 01252 01253 ! ===================== 01254 ! SUBROUTINE HELDSUAREZ 01255 ! ===================== 01256 01257 subroutine heldsuarez 01258 use pumamod 01259 01260 ! Set up the restoration temperature field according to that given 01261 ! in Held & Suarez (1994, Bul. Amer. Meteor. Soc.). 01262 ! Only difference: There is an offset of 1./3. added to the sine of 01263 ! latitude. The reason is that the namelist parameter TGR still 01264 ! represents the global mean of surface restoration temperature. 01265 ! 01266 ! ==> Set TGR=295. in order to get exactly the same restoration 01267 ! temperature field as in Held & Suarez (1994). 01268 ! 01269 ! ==> DTNS in H&S (1994) is the equator-pole difference while 01270 ! PUMA uses DTNS for pole-pole difference. Therefore we use 01271 ! 0.5 * dtns in this subroutine. 01272 01273 01274 ! Produce restoration temperature 01275 01276 real :: zp0,z3 01277 real :: zdtc,zdtv,zsip 01278 01279 zp0 = 100000. 01280 z3 = 1.0 / 3.0 01281 if (nhelsua == 1 .or. nhelsua == 2) then 01282 do jlev=1,nlev 01283 do jlat=1,nlat 01284 zsip = sigma(jlev)*PSURF/zp0 01285 zdtc = dtep * (sid(jlat)**2-z3) + dtzz * log(zsip) * csq(jlat) 01286 zdtv = zdtc - 0.5 * dtns * sid(jlat) 01287 gtc(:,jlat,jlev) = max(ttp, (tgr - zdtc) * zsip**AKAP) 01288 gtv(:,jlat,jlev) = max(ttp, (tgr - zdtv) * zsip**AKAP) 01289 enddo 01290 enddo 01291 gtv(:,:,:) = gtv(:,:,:) - gtc(:,:,:) 01292 endif 01293 01294 ! Produce reciprocal of damping time scale for T [1/sec] 01295 01296 ! tauta = 1.0 / (TWOPI * tauta) 01297 ! tauts = 1.0 / (TWOPI * tauts) 01298 01299 rtauta_dim = 1.0 /(tauta*sid_day) 01300 rtauts_dim = 1.0 /(tauts*sid_day) 01301 if (nhelsua == 2 .or. nhelsua == 3) then 01302 do jlev=1,nlev 01303 if (sigma(jlev) > zsigb) then 01304 do jlat = 1,nlat 01305 gtdamp(1:nlon,jlat,jlev) = rtauta_dim & 01306 & + (rtauts_dim - rtauta_dim) * ((sigma(jlev) - zsigb) / (1. - zsigb)) & 01307 & * (1. - sid(jlat)**2)**2 01308 enddo 01309 else 01310 gtdamp(:,:,jlev) = rtauta_dim 01311 endif 01312 enddo 01313 endif 01314 01315 return 01316 end 01317 01318 ! ===================== 01319 ! SUBROUTINE RESOLUTION 01320 ! ===================== 01321 01322 subroutine resolution 01323 use pumamod 01324 logical :: lex 01325 namelist /res/ nlat, nlev 01326 01327 nlat = 32 01328 nlev = 10 01329 01330 inquire(file=resolution_namelist,exist=lex) 01331 if (.not. lex) then 01332 resolution_namelist = trim(resolution_namelist) // "_00" 01333 inquire(file=resolution_namelist,exist=lex) 01334 endif 01335 01336 if (lex) then 01337 open(14,file=resolution_namelist) 01338 read(14,res) 01339 close(14) 01340 endif 01341 01342 nlem = nlev - 1 01343 nlep = nlev + 1 01344 nlsq = nlev * nlev 01345 01346 nlon = nlat + nlat ! Longitudes 01347 nlah = nlat / 2 01348 nlpp = nlat 01349 nhpp = nlah 01350 nhor = nlon * nlpp 01351 01352 ntru = (nlon - 1) / 3 01353 ntp1 = ntru + 1 01354 nzom = ntp1 + ntp1 01355 nrsp = (ntru + 1) * (ntru + 2) 01356 ncsp = nrsp / 2 01357 nspp = nrsp 01358 nesp = nspp 01359 01360 return 01361 end 01362 01363 ! ================= 01364 ! SUBROUTINE READNL 01365 ! ================= 01366 01367 subroutine readnl 01368 use pumamod 01369 01370 ! This namelist must be identical to namelist inp in puma.f90 01371 01372 namelist /inp/ & 01373 alpha , alrpv , alrs , disp & 01374 , dorox , doroxs , doroy , doroys , dt , dtep & 01375 , dtns , dtrop , dttrp , dtzz , dvdiff , edgepv & 01376 , epsync & 01377 , flsp0 , flsdp , flsamp , flsoff , horo , kick & 01378 , lat1oro , lat1tgr , lat2oro , lat2tgr & 01379 , lon1oro , lon1tgr , lon2oro , lon2tgr & 01380 , mpstep & 01381 , nafter , ncoeff , ncorrect, ncu , ndel , ndiag & 01382 , ndl , nextout , nfls , ngui , nguidbg & 01383 , nhelsua, nsync & 01384 , ntestgp , nkits & 01385 , nlevt , nmonths , noro , norox , noutput , noutsrv & 01386 , npackgp , npacksp , nreverse& 01387 , nruido , nrun , nselect , nspecsel, nsponge , nsrv & 01388 , nstep , nstop , nstrato , ntgr , nsym , ntspd & 01389 , nvg , nwpd , nwspini , nyears , nyoden & 01390 , oroano , orofac & 01391 , pac , pmaxpv , pspon & 01392 , radpv , restim , rotspd & 01393 , sigmah , sigmax , sponk & 01394 , t0 , t0k , tac , tauta , tauts & 01395 , tdiss , tfrc , tgr , tgrano , ttp & 01396 , nenergy , nentropy, ndheat 01397 01398 nselect(:) = 1 01399 nspecsel(:) = 1 01400 ndl(:) = 0 01401 restim(:) = 15.0 01402 sigmah(:) = 0.0 01403 tfrc(1:nlev) = (/ (0.0,i=1,nlem), 1.0 /) 01404 t0k(:) = 250.0 01405 t0(:) = 250.0 01406 dt(:) = 0.0 01407 01408 open(13,file='ppp_namelist') 01409 read (13,inp) 01410 01411 ! Use predefined Yoden profile ? 01412 01413 if (nlev == 20 .and. nyoden > 0) then 01414 ! noro = 0 ! Don't read orography 01415 ! norox = 2 ! Make idealized orography 01416 ! horo = 500.0 ! Height of idealized orography 01417 alrs = -0.001 ! Stratospheric lapse rate 01418 nreverse = 0 ! No T-gradient reversal at tropopause 01419 ncorrect = 0 ! No T-correction due to orography 01420 idim = min(20,nlev) 01421 if (nyoden == 1) then 01422 t0(1:idim) = t0yod1(1:idim) 01423 dt(1:idim) = dtyod1(1:idim) 01424 else if (nyoden == 3) then 01425 t0(1:idim) = t0yod3(1:idim) 01426 dt(1:idim) = dtyod3(1:idim) 01427 else if (nyoden == 5) then 01428 t0(1:idim) = t0yod5(1:idim) 01429 dt(1:idim) = dtyod5(1:idim) 01430 else if (nyoden == 7) then 01431 t0(1:idim) = t0yod7(1:idim) 01432 dt(1:idim) = dtyod7(1:idim) 01433 else if (nyoden == 8) then 01434 t0(1:idim) = t0yod8(1:idim) 01435 dt(1:idim) = dtyod8(1:idim) 01436 else if (nyoden == 9) then 01437 t0(1:idim) = t0yod9(1:idim) 01438 dt(1:idim) = dtyod9(1:idim) 01439 endif 01440 endif 01441 01442 close(13) 01443 write(*,inp) 01444 01445 return 01446 end 01447 01448 ! ===================== 01449 ! * SET VERTICAL GRID * 01450 ! ===================== 01451 01452 subroutine set_vertical_grid 01453 01454 use pumamod 01455 01456 if (sigmah(nlev) /= 0.0) return ! Already read in from namelist INP 01457 01458 if (nvg == 1) then ! Scinocca & Haynes sigma levels 01459 01460 if (nlevt >= nlev) then ! Security check for 'nlevt' 01461 write (*,*) '*** ERROR *** nlevt >= nlev' 01462 write (*,*) 'Number of levels (nlev): ',nlev 01463 write (*,*) 'Number of tropospheric levels (nlevt): ',nlevt 01464 endif 01465 01466 ! troposphere: linear spacing in sigma 01467 ! stratosphere: linear spacing in log(sigma) 01468 ! after (see their Appendix): 01469 ! Scinocca, J. F. and P. H. Haynes (1998): Dynamical forcing of 01470 ! stratospheric planetary waves by tropospheric baroclinic eddies. 01471 ! J. Atmos. Sci., 55 (14), 2361-2392 01472 01473 ! Here, zsigtran is set to sigma at dtrop (tropopause height for 01474 ! construction of restoration temperature field). If tgr=288.15K, 01475 ! ALR=0.0065K/km and dtrop=11.km, then zsigtran=0.223 (=0.1 in 01476 ! Scinocca and Haynes (1998)). 01477 ! A smoothing of the transition between linear and logarithmic 01478 ! spacing, as noted in Scinocca and Haynes (1998), is not yet 01479 ! implemented. 01480 01481 zsigtran = (1. - ALR * dtrop / tgr)**(GA/(GASCON*ALR)) 01482 zsigmin = 1. - (1. - zsigtran) / real(nlevt) 01483 01484 do jlev=1,nlev 01485 if (jlev == 1) then 01486 sigmah(jlev) = 0.000001 01487 sigmah(jlev) = sigmax 01488 elseif (jlev > 1 .and. jlev < nlev - nlevt) then 01489 sigmah(jlev) = exp((log(sigmax) - log(zsigtran)) & 01490 & / real(nlev - nlevt - 1) * real(nlev - nlevt - jlev) 01491 + log(zsigtran)) 01492 elseif (jlev >= nlev - nlevt .and. jlev < nlev - 1) then 01493 sigmah(jlev) = (zsigtran - zsigmin) / real(nlevt - 1) 01494 * real(nlev - 1 - jlev) + zsigmin 01495 elseif (jlev == nlev - 1) then 01496 sigmah(jlev) = zsigmin 01497 elseif (jlev == nlev) then 01498 sigmah(jlev) = 1. 01499 endif 01500 enddo 01501 return 01502 endif ! (nvg == 1) 01503 01504 if (nvg == 2) then ! Polvani & Kushner sigma levels 01505 inl = int(real(nlev)/(1.0 - sigmax**(1.0/5.0))) 01506 do jlev=1,nlev 01507 sigmah(jlev) = (real(jlev + inl - nlev) / real(inl))**5 01508 enddo 01509 return 01510 endif ! (nvg == 2) 01511 01512 ! Default: equidistant sigma levels 01513 01514 if (nvg == 0) then 01515 do jlev = 1 , nlev 01516 sigmah(jlev) = real(jlev) / real(nlev) 01517 enddo 01518 endif ! (nvg == 0) 01519 01520 return 01521 end 01522 01523 ! ================= 01524 ! SUBROUTINE INITPM 01525 ! ================= 01526 01527 subroutine initpm 01528 use pumamod 01529 01530 real (kind=8) radea,zakk,zzakk 01531 radea = PLARAD_EARTH ! Planet radius in high precision 01532 plavor = EZ * rotspd ! Planetary vorticity 01533 01534 ! ************************************************************* 01535 ! * carries out all initialisation of model prior to running. * 01536 ! * major sections identified with comments. * 01537 ! * this s/r sets the model parameters and all resolution * 01538 ! * dependent quantities. * 01539 ! ************************************************************* 01540 01541 ! ********************* 01542 ! * set vertical grid * 01543 ! ********************* 01544 01545 call set_vertical_grid 01546 01547 dsigma(1 ) = sigmah(1) 01548 dsigma(2:nlev) = sigmah(2:nlev) - sigmah(1:NLEM) 01549 01550 rdsig(:) = 0.5 / dsigma(:) 01551 01552 sigma(1 ) = 0.5 * sigmah(1) 01553 sigma(2:nlev) = 0.5 * (sigmah(1:NLEM) + sigmah(2:nlev)) 01554 01555 ! annual cycle period and phase in timesteps 01556 01557 if (tac > 0.0) tac = TWOPI / (ntspd * tac) 01558 pac = pac * ntspd 01559 01560 ! compute internal diffusion parameter 01561 01562 jdelh = ndel/2 01563 if (tdiss > 0.0) then 01564 zakk = WW*(radea**ndel)/(TWOPI*tdiss*((ntru*(ntru+1.))**jdelh)) 01565 else 01566 zakk = 0.0 01567 endif 01568 zzakk = zakk / (WW*(radea**ndel)) 01569 01570 ! set coefficients which depend on wavenumber 01571 01572 zrsq2 = 1.0 / sqrt(2.0) 01573 01574 jr=-1 01575 do jm=0,ntru 01576 do jn=jm,ntru 01577 jr=jr+2 01578 ji=jr+1 01579 spnorm(jr)=zrsq2 01580 spnorm(ji)=zrsq2 01581 enddo 01582 zrsq2=-zrsq2 01583 enddo 01584 01585 return 01586 end 01587 01588 subroutine printparameter 01589 use pumamod 01590 01591 print 8000 01592 print 8050 01593 print 8000 01594 print 8010,nlev 01595 print 8020,ntru 01596 print 8030,nlat 01597 print 8040,nlon 01598 print 8000 01599 print 8120 01600 return 01601 8000 format(' *****************************************') 01602 8010 format(' * nlev = ',i6,' Number of levels *') 01603 8020 format(' * ntru = ',i6,' Triangular truncation *') 01604 8030 format(' * nlat = ',i6,' Number of latitudes *') 01605 8040 format(' * nlon = ',i6,' Number of longitues *') 01606 8050 format(' * PPP - Puma Pre Processor *') 01607 8120 format(/) 01608 end 01609 01610 01611 ! =================== 01612 ! SUBROUTINE GPROFILE 01613 ! =================== 01614 subroutine gprofile(ptgr,prgrad,pgpm) 01615 use pumamod 01616 01617 ! ************************************************************* 01618 ! * Set up the restoration temperature profiles for gradient * 01619 ! * modes DTNS - mode[0,1] and DTEP - mode[0,2] * 01620 ! * The lapse rate of ALR K/m is assumed under the tropopause * 01621 ! * and zero above. The tropopause is defined by <dtrop>. * 01622 ! * The profile is a sine wave with 0 at tropopause sigma and * 01623 ! * 1 at sigma = 1. * 01624 ! ************************************************************* 01625 01626 dimension prgrad(nlev) 01627 01628 ztpheight = dtrop - pgpm ! Tropopause height over ground 01629 ztptemp = tgr - ALR * dtrop ! Tropopause temperature 01630 ztps = (ztptemp/ptgr)**(GA/(ALR*GASCON)) ! Tropoause sigma 01631 01632 do jlev = 1 , nlev 01633 prgrad(jlev) = sin(0.5*PI*(sigma(jlev)-ztps)/(1.0-ztps)) 01634 if (sigma(jlev) < ztps) prgrad(jlev) = 0.0 01635 enddo 01636 01637 return 01638 end 01639 01640 ! =================== 01641 ! SUBROUTINE TPROFILE 01642 ! =================== 01643 subroutine tprofile(ptgr,prof,pgpm) 01644 use pumamod 01645 01646 ! ************************************************************* 01647 ! * Set up the restoration temperature profile for one column * 01648 ! * The temperature at sigma = 1 is <ptgr>, entered in kelvin * 01649 ! * The lapse rate of ALR K/m is assumed under the tropopause * 01650 ! * and zero above. The tropopause is defined by <ztpheight>. * 01651 ! * The smoothing ot the tropopause depends on <dttrp>. * 01652 ! ************************************************************* 01653 01654 dimension prof(nlev) ! Resulting temperature profile [K] 01655 01656 zsigprev = 1.0 ! sigma value 01657 ztprev = ptgr ! Temperature [K] 01658 zzprev = 0.0 ! Height [m] 01659 ztpheight = dtrop - pgpm ! Tropopause height over ground 01660 ztptemp = tgr - ALR * dtrop ! Tropopause temperature 01661 zalr = (ptgr - ztptemp) / ztpheight 01662 01663 do jlev = nlev , 1 , -1 ! from bottom to top of atmosphere 01664 zlogsig = GASCON / GA * log(zsigprev / sigma(jlev)) 01665 zzp = zzprev + ztprev * zlogsig 01666 ztp=ztptemp+sqrt((.5*zalr*(zzp-ztpheight))**2+dttrp**2) 01667 ztp=ztp-.5*zalr*(zzp-ztpheight) 01668 ztpm=.5*(ztprev+ztp) 01669 01670 zzpp = zzprev + ztpm * zlogsig 01671 ztpp=ztptemp+sqrt((.5*zalr*(zzpp-ztpheight))**2+dttrp**2) 01672 ztpp=ztpp-.5*zalr*(zzpp-ztpheight) 01673 01674 prof(jlev)=ztpp 01675 zzprev=zzprev + 0.5 * (ztpp+ztprev) * zlogsig 01676 ztprev=ztpp 01677 zsigprev=sigma(jlev) 01678 enddo 01679 return 01680 end 01681 01682 ! ====================== 01683 ! SUBROUTINE ppp_write_i 01684 ! ====================== 01685 01686 subroutine ppp_write_i(yvarname,nvals,ivals) 01687 use pumamod 01688 01689 character(*) :: yvarname 01690 integer :: nvals 01691 integer :: ivals(nvals) 01692 01693 if (mypid == NROOT) then 01694 write(95,'("[",A,"]")') trim(yvarname) 01695 write(95,'(I4)') nvals 01696 write(95,'(I6)') ivals(:) 01697 endif 01698 return 01699 end 01700 01701 ! ====================== 01702 ! SUBROUTINE ppp_write_r 01703 ! ====================== 01704 01705 subroutine ppp_write_r(yvarname,nvals,pvals) 01706 use pumamod 01707 01708 character(*) :: yvarname 01709 integer :: nvals 01710 real :: pvals(nvals) 01711 01712 if (mypid == NROOT) then 01713 write(95,'("[",A,"]")') trim(yvarname) 01714 write(95,'(I4)') nvals 01715 write(95,'(E14.8)') pvals(:) 01716 endif 01717 return 01718 end 01719 01720 ! ================ 01721 ! SUBROUTINE yoden 01722 ! ================ 01723 01724 subroutine yoden 01725 use pumamod 01726 ! 01727 do jlev=1,nlev 01728 do jlat=1,nlat/2 01729 zlat=dla(jlat) 01730 ztres=t0(jlev)+dt(jlev)/2.*(cos(2.*zlat)-1./3.) 01731 do jlon=1,nlon 01732 gtc(jlon,jlat,jlev)=ztres 01733 enddo 01734 enddo 01735 enddo 01736 ! 01737 do jlat=1,nlat/2 01738 j2=nlat+1-jlat 01739 gtc(:,j2,:)=gtc(:,jlat,:) 01740 enddo 01741 01742 write (*,*) ' Computed Yoden profile',nyoden 01743 ! 01744 return 01745 end 01746 ! 01747 ! ======================= 01748 ! SUBROUTINE PRINTPROFILE 01749 ! ======================= 01750 01751 subroutine printprofile 01752 use pumamod 01753 01754 ! ********************************** 01755 ! * write out vertical information * 01756 ! ********************************** 01757 01758 dimension ztr(nlev+1) 01759 01760 ztr(nlev+1) = tgr 01761 01762 write(*,9001) 01763 write(*,9002) 01764 write(*,9003) 01765 write(*,9002) 01766 01767 do jlev=1,nlev 01768 ztr(jlev) = sum(gtc(:,:,jlev)) / (nlon * nlat) 01769 enddo 01770 01771 do jlev=1,nlev 01772 write(*,9004) jlev,sigma(jlev),ztr(jlev),ztr(jlev+1)-ztr(jlev),gra(jlev) 01773 enddo 01774 01775 write(*,9002) 01776 write(*,9001) 01777 return 01778 9001 format(/) 01779 9002 format(1x,45('*')) 01780 9003 format(' * Lv * Sigma Restor-T Delta-T Vfact *') 01781 9004 format(' *',i3,' * ',4f9.3,' *') 01782 end 01783 01784 ! ===================== 01785 ! * SUBROUTINE LEGPRI * 01786 ! ===================== 01787 01788 subroutine legpri 01789 use pumamod 01790 01791 write (*,231) 01792 write (*,232) 01793 write (*,233) 01794 write (*,232) 01795 do 14 jlat = 1 , nlat 01796 zalat = dla(jlat)*180.0/PI 01797 write (*,234) jlat,zalat,csq(jlat),gwd(jlat) 01798 14 continue 01799 write (*,232) 01800 write (*,231) 01801 return 01802 231 format(/) 01803 232 format(1x,36('*')) 01804 233 format(' * No * Lat * csq weight *') 01805 234 format(' *',i3,' *',f6.1,' *',2f10.4,' *') 01806 end 01807 01808 01809 ! ================= 01810 ! SUBROUTINE INILAT 01811 ! ================= 01812 01813 subroutine inilat 01814 use pumamod 01815 character(1) :: ch 01816 01817 ch = 'N' 01818 do jlat = 1 , nlat 01819 csq(jlat) = 1.0 - sid(jlat) * sid(jlat) 01820 dla(jlat) = asin(sid(jlat)) 01821 enddo 01822 do jlat = 1 , nlat/2 01823 ideg = nint(180.0/PI * asin(sid(jlat))) 01824 write(chlat(jlat),'(i2,a1)') ideg,'N' 01825 write(chlat(nlat+1-jlat),'(i2,a1)') ideg,'S' 01826 enddo 01827 return 01828 end 01829 01830 01831 ! ================= 01832 ! SUBROUTINE SETZT2 01833 ! ================= 01834 subroutine setzt2 01835 use pumamod 01836 01837 ! US standard atmosphere (1976): 01838 parameter(INL = 7) ! number of defined layers 01839 dimension zzus(0:INL) ! height of interfaces between layers 01840 dimension zlus(INL) ! temperature lapse rates of layers 01841 dimension zpus(0:INL) ! pressure at interfaces between layers 01842 dimension ztus(0:INL) ! temperature at interfaces between layers 01843 01844 dimension ztrs(nlev) ! Mean profile 01845 dimension ztpv(nlev) ! Vertical profile of stratospheric polar 01846 ! vortex forcing 01847 dimension zdtep(nlat) 01848 dimension zdtns(nlat) 01849 dimension zff(nlev) 01850 dimension zfw(nlat,nlev) 01851 dimension zfsph(nlat,nlev) 01852 dimension zqc(nlat,nlev) 01853 dimension zphi(nlat) 01854 dimension zgr1(nlon,nlat,nlev) 01855 dimension zgr2(nlon,nlat,nlev) 01856 dimension zfls(nlat,nlev) 01857 01858 real :: zp 01859 real :: zsigtp 01860 real :: pref 01861 real :: zpmaxsph 01862 01863 sr1(:,:) = 0.0 ! NESP,nlev 01864 sr2(:,:) = 0.0 ! NESP,nlev 01865 01866 ! 1. Mean vertical profile (MVP), approx. US standard atmosphere 01867 01868 zzus(0) = 0. 01869 zzus(1) = dtrop ! US standard atmosphere: zzus(1) = 11000. 01870 zzus(2) = 20000. 01871 zzus(3) = 32000. 01872 zzus(4) = 47000. 01873 zzus(5) = 51000. 01874 zzus(6) = 71000. 01875 zzus(7) = 84852. 01876 zlus(1) = ALR ! US standard atmosphere: zlus(1) = 0.0065 01877 zlus(2) = 0.0 01878 zlus(3) = -0.001 01879 zlus(4) = -0.0028 01880 zlus(5) = 0.0 01881 zlus(6) = 0.0028 01882 zlus(7) = 0.002 01883 01884 ! calculation of pressure and temperature at layer interfaces 01885 01886 zpus(0) = PSURF ! US standard atmosphere: zpus(0) = 1013.25 hPa 01887 ztus(0) = tgr ! US standard atmosphere: ztus(0) = 288.15 K 01888 01889 do ji=1,INL 01890 ztus(ji) = ztus(ji-1) - zlus(ji) * (zzus(ji) - zzus(ji-1)) 01891 if (zlus(ji) == 0.) then 01892 zpus(ji) = zpus(ji-1) * exp(-GA * (zzus(ji) - zzus(ji-1)) & 01893 & / (GASCON * ztus(ji-1))) 01894 else 01895 zpus(ji) = zpus(ji-1) & 01896 & * (ztus(ji) / ztus(ji-1))**(GA/(GASCON*zlus(ji))) 01897 endif 01898 enddo 01899 01900 ! calculation of temperature on given sigma full levels, sigma(1:nlev) 01901 do jlev=nlev,1,-1 01902 zp = sigma(jlev)*PSURF 01903 if (zp <= zpus(0) .and. zp > zpus(1)) then 01904 ztrs(jlev) = ztus(0) * (zp / zpus(0))**(GASCON*zlus(1)/GA) 01905 elseif (zp <= zpus(1) .and. zp > zpus(2)) then 01906 ztrs(jlev) = ztus(1) * (zp / zpus(1))**(GASCON*zlus(2)/GA) 01907 elseif (zp <= zpus(2) .and. zp > zpus(3)) then 01908 ztrs(jlev) = ztus(2) * (zp / zpus(2))**(GASCON*zlus(3)/GA) 01909 elseif (zp <= zpus(3) .and. zp > zpus(4)) then 01910 ztrs(jlev) = ztus(3) * (zp / zpus(3))**(GASCON*zlus(4)/GA) 01911 elseif (zp <= zpus(4) .and. zp > zpus(5)) then 01912 ztrs(jlev) = ztus(4) * (zp / zpus(4))**(GASCON*zlus(5)/GA) 01913 elseif (zp <= zpus(5) .and. zp > zpus(6)) then 01914 ztrs(jlev) = ztus(5) * (zp / zpus(5))**(GASCON*zlus(6)/GA) 01915 elseif (zp <= zpus(6) .and. zp > zpus(7)) then 01916 ztrs(jlev) = ztus(6) * (zp / zpus(6))**(GASCON*zlus(7)/GA) 01917 else 01918 ztrs(jlev) = ztus(7) 01919 endif 01920 enddo 01921 01922 ! 2. Symmetric equator-pole forcing mode (DTEP) and 01923 ! 3. Asymmetric Npole-Spole forcing mode (DTNS) 01924 01925 ! sid(nlat) is sine of latitude, taking into account the nonequally 01926 ! spaced Gaussian latitudes. 01927 do jlat=1,nlat 01928 zdtep(jlat) = -dtep * (sid(jlat)**2 - 1./3.) 01929 zdtns(jlat) = dtns * sid(jlat) / 2. 01930 enddo 01931 01932 ! 4. Factor modulating the DTEP and DTNS modes (f) 01933 01934 zsigtp = zpus(1)/zpus(0) ! sigma at tropopause 01935 zff(:) = 0. 01936 do jlev=1,nlev 01937 if (sigma(jlev) > zsigtp) then 01938 zff(jlev) = sin(0.5*PI * (sigma(jlev) - zsigtp) & 01939 & / (1. - zsigtp)) 01940 endif 01941 enddo 01942 01943 ! 5. Vertical profile of stratospheric polar vortex forcing 01944 01945 ztpv(:) = 0. 01946 if (pmaxpv == 0.) then 01947 do jlev=1,nlev 01948 if (sigma(jlev) <= zsigtp) then 01949 ztpv(jlev) = ztus(1) * (sigma(jlev)*PSURF / zpus(1)) & 01950 & **(GASCON*alrpv/GA) 01951 endif 01952 enddo 01953 elseif (pmaxpv > 0.) then 01954 do jlev=1,nlev 01955 if (sigma(jlev) <= pmaxpv/PSURF) then 01956 ztpv(jlev) = ztus(1) * (sigma(jlev)*PSURF / pmaxpv) & 01957 & **(GASCON*alrpv/GA) 01958 else 01959 ztpv(jlev) = ztus(1) 01960 endif 01961 enddo 01962 endif 01963 01964 ! 6. Factor confining the stratosph. polar vortex to high latitudes 01965 01966 zphi(:) = dla(:) * 180. / PI 01967 zfw(:,:) = 0. 01968 if (edgepv > 0.) then 01969 do jlev=1,nlev 01970 if (sigma(jlev) <= pmaxpv/PSURF) then 01971 do jlat=1,nlat 01972 zfw(jlat,jlev) = & 01973 & 0.5 * (1. - tanh((radpv - zphi(jlat)) / edgepv)) 01974 enddo 01975 endif 01976 enddo 01977 endif 01978 01979 ! 7. Lower stratospheric forcing 01980 01981 zfls(:,:) = 0. 01982 if (nfls == 1) then 01983 do jlev=1,nlev 01984 zp =sigma(jlev) * PSURF 01985 do jlat=1,nlat 01986 if (zp > flsp0-flsdp .and. zp < flsp0+flsdp) then 01987 zfls(jlat,jlev) = cos(0.5 * PI * (zp - flsp0) / flsdp)& 01988 & * (flsoff - flsamp * cos(2. * zphi(jlat) * PI / 180.)) 01989 endif 01990 enddo 01991 enddo 01992 endif 01993 01994 ! construct restoration temperature field 01995 01996 do jlev=1,nlev 01997 do jlat=1,nlat 01998 zgr1(:,jlat,jlev) = ((1. - zfw(jlat,jlev)) * ztrs(jlev) & 01999 & + zfw(jlat,jlev) * ztpv(jlev) + zff(jlev) * zdtep(jlat) & 02000 & + (1. - zfw(jlat,jlev)) * zfls(jlat,jlev) - t0k(jlev)) / CT 02001 zgr2(:,jlat,jlev) = (zff(jlev) * zdtns(jlat)) / CT 02002 enddo 02003 enddo 02004 02005 do jlev = 1 , nlev 02006 gtc(:,:,jlev) = t0k(jlev) + CT * zgr1(:,:,jlev) 02007 gtv(:,:,jlev) = CT * zgr2(:,:,jlev) 02008 enddo 02009 02010 ! ---------- test output to control T_r field ---------- 02011 open(112,file='tr_test.srv',form='unformatted') 02012 do jlev=1,nlev 02013 ip = int(sigma(jlev) * PSURF * 1000.) 02014 write(112) 121, ip, 0000, 00, nlon, nlat, 0, 0 02015 write(112) t0k(jlev) + (zgr1(:,:,jlev) + zgr2(:,:,jlev)) * CT 02016 enddo 02017 close(112) 02018 ! ---------- test output to control T_r field ---------- 02019 02020 print *,'**************************************************' 02021 print *,'* Restoration Temperature set up for aqua planet *' 02022 print *,'* including stratosphere and polar vortex *' 02023 print *,'**************************************************' 02024 return 02025 end
