[234] | 1 | SUBROUTINE physiq( |
---|
[226] | 2 | $ ngrid,nlayer,nq |
---|
| 3 | $ ,firstcall,lastcall |
---|
| 4 | $ ,pday,ptime,ptimestep |
---|
| 5 | $ ,pplev,pplay,pphi |
---|
| 6 | $ ,pu,pv,pt,pq |
---|
| 7 | $ ,pw |
---|
| 8 | $ ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn |
---|
[86] | 9 | #ifdef MESOSCALE |
---|
[226] | 10 | #include "meso_inc/meso_inc_invar.F" |
---|
[86] | 11 | #endif |
---|
[226] | 12 | $ ) |
---|
[42] | 13 | |
---|
[544] | 14 | |
---|
[42] | 15 | IMPLICIT NONE |
---|
| 16 | c======================================================================= |
---|
| 17 | c |
---|
| 18 | c subject: |
---|
| 19 | c -------- |
---|
| 20 | c |
---|
| 21 | c Organisation of the physical parametrisations of the LMD |
---|
| 22 | c martian atmospheric general circulation model. |
---|
| 23 | c |
---|
| 24 | c The GCM can be run without or with tracer transport |
---|
| 25 | c depending on the value of Logical "tracer" in file "callphys.def" |
---|
| 26 | c Tracers may be water vapor, ice OR chemical species OR dust particles |
---|
| 27 | c |
---|
| 28 | c SEE comments in initracer.F about numbering of tracer species... |
---|
| 29 | c |
---|
| 30 | c It includes: |
---|
| 31 | c |
---|
| 32 | c 1. Initialization: |
---|
| 33 | c 1.1 First call initializations |
---|
| 34 | c 1.2 Initialization for every call to physiq |
---|
| 35 | c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff. |
---|
| 36 | c 2. Compute radiative transfer tendencies |
---|
| 37 | c (longwave and shortwave) for CO2 and aerosols. |
---|
| 38 | c 3. Gravity wave and subgrid scale topography drag : |
---|
| 39 | c 4. Vertical diffusion (turbulent mixing): |
---|
| 40 | c 5. Convective adjustment |
---|
| 41 | c 6. Condensation and sublimation of carbon dioxide. |
---|
| 42 | c 7. TRACERS : |
---|
| 43 | c 7a. water and water ice |
---|
| 44 | c 7b. call for photochemistry when tracers are chemical species |
---|
| 45 | c 7c. other scheme for tracer (dust) transport (lifting, sedimentation) |
---|
| 46 | c 7d. updates (CO2 pressure variations, surface budget) |
---|
| 47 | c 8. Contribution to tendencies due to thermosphere |
---|
| 48 | c 9. Surface and sub-surface temperature calculations |
---|
| 49 | c 10. Write outputs : |
---|
| 50 | c - "startfi", "histfi" (if it's time) |
---|
| 51 | c - Saving statistics (if "callstats = .true.") |
---|
| 52 | c - Dumping eof (if "calleofdump = .true.") |
---|
| 53 | c - Output any needed variables in "diagfi" |
---|
| 54 | c 11. Diagnostic: mass conservation of tracers |
---|
| 55 | c |
---|
| 56 | c author: |
---|
| 57 | c ------- |
---|
| 58 | c Frederic Hourdin 15/10/93 |
---|
| 59 | c Francois Forget 1994 |
---|
| 60 | c Christophe Hourdin 02/1997 |
---|
| 61 | c Subroutine completly rewritten by F.Forget (01/2000) |
---|
| 62 | c Introduction of the photochemical module: S. Lebonnois (11/2002) |
---|
| 63 | c Introduction of the thermosphere module: M. Angelats i Coll (2002) |
---|
| 64 | c Water ice clouds: Franck Montmessin (update 06/2003) |
---|
| 65 | c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009) |
---|
| 66 | c Nb: See callradite.F for more information. |
---|
[234] | 67 | c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags |
---|
[414] | 68 | c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization |
---|
[42] | 69 | c |
---|
| 70 | c arguments: |
---|
| 71 | c ---------- |
---|
| 72 | c |
---|
| 73 | c input: |
---|
| 74 | c ------ |
---|
| 75 | c ecri period (in dynamical timestep) to write output |
---|
| 76 | c ngrid Size of the horizontal grid. |
---|
| 77 | c All internal loops are performed on that grid. |
---|
| 78 | c nlayer Number of vertical layers. |
---|
| 79 | c nq Number of advected fields |
---|
| 80 | c firstcall True at the first call |
---|
| 81 | c lastcall True at the last call |
---|
| 82 | c pday Number of days counted from the North. Spring |
---|
| 83 | c equinoxe. |
---|
| 84 | c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT |
---|
| 85 | c ptimestep timestep (s) |
---|
| 86 | c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa) |
---|
| 87 | c pplev(ngrid,nlayer+1) intermediate pressure levels (pa) |
---|
| 88 | c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2) |
---|
| 89 | c pu(ngrid,nlayer) u component of the wind (ms-1) |
---|
| 90 | c pv(ngrid,nlayer) v component of the wind (ms-1) |
---|
| 91 | c pt(ngrid,nlayer) Temperature (K) |
---|
| 92 | c pq(ngrid,nlayer,nq) Advected fields |
---|
[330] | 93 | c pudyn(ngrid,nlayer) | |
---|
| 94 | c pvdyn(ngrid,nlayer) | Dynamical temporal derivative for the |
---|
| 95 | c ptdyn(ngrid,nlayer) | corresponding variables |
---|
| 96 | c pqdyn(ngrid,nlayer,nq) | |
---|
[42] | 97 | c pw(ngrid,?) vertical velocity |
---|
| 98 | c |
---|
| 99 | c output: |
---|
| 100 | c ------- |
---|
| 101 | c |
---|
[330] | 102 | c pdu(ngrid,nlayermx) | |
---|
| 103 | c pdv(ngrid,nlayermx) | Temporal derivative of the corresponding |
---|
| 104 | c pdt(ngrid,nlayermx) | variables due to physical processes. |
---|
| 105 | c pdq(ngrid,nlayermx,nqmx) | |
---|
| 106 | c pdpsrf(ngrid) | |
---|
[42] | 107 | c tracerdyn call tracer in dynamical part of GCM ? |
---|
| 108 | |
---|
| 109 | c |
---|
| 110 | c======================================================================= |
---|
| 111 | c |
---|
| 112 | c 0. Declarations : |
---|
| 113 | c ------------------ |
---|
| 114 | |
---|
| 115 | #include "dimensions.h" |
---|
| 116 | #include "dimphys.h" |
---|
| 117 | #include "comgeomfi.h" |
---|
| 118 | #include "surfdat.h" |
---|
[226] | 119 | #include "comsoil.h" |
---|
[42] | 120 | #include "comdiurn.h" |
---|
| 121 | #include "callkeys.h" |
---|
| 122 | #include "comcstfi.h" |
---|
| 123 | #include "planete.h" |
---|
| 124 | #include "comsaison.h" |
---|
| 125 | #include "control.h" |
---|
| 126 | #include "dimradmars.h" |
---|
| 127 | #include "comg1d.h" |
---|
| 128 | #include "tracer.h" |
---|
| 129 | #include "nlteparams.h" |
---|
| 130 | |
---|
| 131 | #include "chimiedata.h" |
---|
| 132 | #include "param.h" |
---|
| 133 | #include "param_v3.h" |
---|
| 134 | #include "conc.h" |
---|
| 135 | |
---|
| 136 | #include "netcdf.inc" |
---|
| 137 | |
---|
[234] | 138 | #include "slope.h" |
---|
| 139 | |
---|
[86] | 140 | #ifdef MESOSCALE |
---|
[42] | 141 | #include "wrf_output_2d.h" |
---|
| 142 | #include "wrf_output_3d.h" |
---|
[226] | 143 | #include "advtrac.h" !!! this is necessary for tracers (in dyn3d) |
---|
| 144 | #include "meso_inc/meso_inc_var.F" |
---|
[86] | 145 | #endif |
---|
[42] | 146 | |
---|
| 147 | c Arguments : |
---|
| 148 | c ----------- |
---|
| 149 | |
---|
| 150 | c inputs: |
---|
| 151 | c ------- |
---|
| 152 | INTEGER ngrid,nlayer,nq |
---|
| 153 | REAL ptimestep |
---|
| 154 | REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer) |
---|
| 155 | REAL pphi(ngridmx,nlayer) |
---|
| 156 | REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer) |
---|
| 157 | REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq) |
---|
| 158 | REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d |
---|
| 159 | REAL zh(ngridmx,nlayermx) ! potential temperature (K) |
---|
| 160 | LOGICAL firstcall,lastcall |
---|
[226] | 161 | |
---|
[42] | 162 | REAL pday |
---|
| 163 | REAL ptime |
---|
| 164 | logical tracerdyn |
---|
| 165 | |
---|
| 166 | c outputs: |
---|
| 167 | c -------- |
---|
| 168 | c physical tendencies |
---|
| 169 | REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer) |
---|
| 170 | REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq) |
---|
| 171 | REAL pdpsrf(ngridmx) ! surface pressure tendency |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | c Local saved variables: |
---|
| 175 | c ---------------------- |
---|
| 176 | c aerosol (dust or ice) extinction optical depth at reference wavelength |
---|
| 177 | c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of |
---|
| 178 | c aerosol optical properties : |
---|
| 179 | REAL aerosol(ngridmx,nlayermx,naerkind) |
---|
| 180 | |
---|
| 181 | INTEGER day_ini ! Initial date of the run (sol since Ls=0) |
---|
| 182 | INTEGER icount ! counter of calls to physiq during the run. |
---|
| 183 | REAL tsurf(ngridmx) ! Surface temperature (K) |
---|
| 184 | REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K) |
---|
| 185 | REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2) |
---|
| 186 | REAL albedo(ngridmx,2) ! Surface albedo in each solar band |
---|
| 187 | REAL emis(ngridmx) ! Thermal IR surface emissivity |
---|
| 188 | REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1) |
---|
| 189 | REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2) |
---|
| 190 | REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2) |
---|
| 191 | REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1) |
---|
| 192 | REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2) |
---|
| 193 | REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2) |
---|
[283] | 194 | REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy |
---|
| 195 | |
---|
| 196 | REAL watercapflag(ngridmx) ! water cap flag |
---|
[42] | 197 | |
---|
| 198 | c Variables used by the water ice microphysical scheme: |
---|
| 199 | REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m) |
---|
| 200 | REAL nuice(ngridmx,nlayermx) ! Estimated effective variance |
---|
| 201 | ! of the size distribution |
---|
[358] | 202 | real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius |
---|
| 203 | real rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3) |
---|
[520] | 204 | REAL surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3, if photochemistry) |
---|
| 205 | REAL surfice(ngridmx,nlayermx) ! ice surface area (micron2/cm3, if photochemistry) |
---|
[42] | 206 | |
---|
[234] | 207 | c Variables used by the slope model |
---|
| 208 | REAL sl_ls, sl_lct, sl_lat |
---|
| 209 | REAL sl_tau, sl_alb, sl_the, sl_psi |
---|
| 210 | REAL sl_fl0, sl_flu |
---|
| 211 | REAL sl_ra, sl_di0 |
---|
| 212 | REAL sky |
---|
| 213 | |
---|
[42] | 214 | SAVE day_ini, icount |
---|
| 215 | SAVE aerosol, tsurf,tsoil |
---|
| 216 | SAVE co2ice,albedo,emis, q2 |
---|
| 217 | SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf |
---|
| 218 | |
---|
| 219 | REAL stephan |
---|
| 220 | DATA stephan/5.67e-08/ ! Stephan Boltzman constant |
---|
| 221 | SAVE stephan |
---|
| 222 | |
---|
| 223 | c Local variables : |
---|
| 224 | c ----------------- |
---|
| 225 | |
---|
| 226 | REAL CBRT |
---|
| 227 | EXTERNAL CBRT |
---|
| 228 | |
---|
| 229 | CHARACTER*80 fichier |
---|
| 230 | INTEGER l,ig,ierr,igout,iq,i, tapphys |
---|
| 231 | |
---|
| 232 | REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2) |
---|
| 233 | REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2) |
---|
| 234 | REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2) |
---|
| 235 | REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2) |
---|
| 236 | REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa |
---|
| 237 | REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point |
---|
| 238 | REAL zls ! solar longitude (rad) |
---|
| 239 | REAL zday ! date (time since Ls=0, in martian days) |
---|
| 240 | REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers |
---|
| 241 | REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries |
---|
| 242 | REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic) |
---|
| 243 | |
---|
| 244 | c Tendancies due to various processes: |
---|
| 245 | REAL dqsurf(ngridmx,nqmx) |
---|
| 246 | REAL zdtlw(ngridmx,nlayermx) ! (K/s) |
---|
| 247 | REAL zdtsw(ngridmx,nlayermx) ! (K/s) |
---|
| 248 | REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area |
---|
| 249 | REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area |
---|
| 250 | REAL zdtnirco2(ngridmx,nlayermx) ! (K/s) |
---|
| 251 | REAL zdtnlte(ngridmx,nlayermx) ! (K/s) |
---|
| 252 | REAL zdtsurf(ngridmx) ! (K/s) |
---|
| 253 | REAL zdtcloud(ngridmx,nlayermx) |
---|
| 254 | REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2) |
---|
| 255 | REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s) |
---|
| 256 | REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2) |
---|
| 257 | REAL zdhadj(ngridmx,nlayermx) ! (K/s) |
---|
| 258 | REAL zdtgw(ngridmx,nlayermx) ! (K/s) |
---|
| 259 | REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2) |
---|
| 260 | REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx) |
---|
| 261 | REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx) |
---|
| 262 | |
---|
| 263 | REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx) |
---|
| 264 | REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx) |
---|
| 265 | REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx) |
---|
| 266 | REAL zdqadj(ngridmx,nlayermx,nqmx) |
---|
| 267 | REAL zdqc(ngridmx,nlayermx,nqmx) |
---|
| 268 | REAL zdqcloud(ngridmx,nlayermx,nqmx) |
---|
| 269 | REAL zdqscloud(ngridmx,nqmx) |
---|
| 270 | REAL zdqchim(ngridmx,nlayermx,nqmx) |
---|
| 271 | REAL zdqschim(ngridmx,nqmx) |
---|
| 272 | |
---|
| 273 | REAL zdteuv(ngridmx,nlayermx) ! (K/s) |
---|
| 274 | REAL zdtconduc(ngridmx,nlayermx) ! (K/s) |
---|
| 275 | REAL zdumolvis(ngridmx,nlayermx) |
---|
| 276 | REAL zdvmolvis(ngridmx,nlayermx) |
---|
| 277 | real zdqmoldiff(ngridmx,nlayermx,nqmx) |
---|
| 278 | |
---|
| 279 | c Local variable for local intermediate calcul: |
---|
| 280 | REAL zflubid(ngridmx) |
---|
| 281 | REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx) |
---|
| 282 | REAL zdum1(ngridmx,nlayermx) |
---|
| 283 | REAL zdum2(ngridmx,nlayermx) |
---|
| 284 | REAL ztim1,ztim2,ztim3, z1,z2 |
---|
| 285 | REAL ztime_fin |
---|
| 286 | REAL zdh(ngridmx,nlayermx) |
---|
| 287 | INTEGER length |
---|
| 288 | PARAMETER (length=100) |
---|
| 289 | |
---|
| 290 | c local variables only used for diagnostic (output in file "diagfi" or "stats") |
---|
| 291 | c ----------------------------------------------------------------------------- |
---|
| 292 | REAL ps(ngridmx), zt(ngridmx,nlayermx) |
---|
| 293 | REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx) |
---|
| 294 | REAL zq(ngridmx,nlayermx,nqmx) |
---|
| 295 | REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx) |
---|
| 296 | character*2 str2 |
---|
| 297 | character*5 str5 |
---|
| 298 | real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx) |
---|
[358] | 299 | REAL tauscaling(ngridmx) ! Convertion factor for qdust and Ndust |
---|
| 300 | SAVE tauscaling ! in case iradia NE 1 |
---|
[42] | 301 | real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m) |
---|
| 302 | integer igmin, lmin |
---|
| 303 | logical tdiag |
---|
| 304 | |
---|
| 305 | real co2col(ngridmx) ! CO2 column |
---|
| 306 | REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx) |
---|
| 307 | REAL zstress(ngrid), cd |
---|
| 308 | real hco2(nqmx),tmean, zlocal(nlayermx) |
---|
| 309 | real rho(ngridmx,nlayermx) ! density |
---|
| 310 | real vmr(ngridmx,nlayermx) ! volume mixing ratio |
---|
[520] | 311 | !real colden(ngridmx,nqmx) ! vertical column !FL |
---|
[42] | 312 | REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2) |
---|
| 313 | REAL icetot(ngridmx) ! Total mass of water ice (kg/m2) |
---|
[520] | 314 | REAL ccntot(ngridmx) ! Total number of ccn (nbr/m2) |
---|
[42] | 315 | REAL rave(ngridmx) ! Mean water ice effective radius (m) |
---|
| 316 | REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1 |
---|
| 317 | REAL tauTES(ngridmx) ! column optical depth at 825 cm-1 |
---|
| 318 | REAL Qabsice ! Water ice absorption coefficient |
---|
[520] | 319 | REAL taucloudtes(ngridmx)! Cloud opacity at infrared |
---|
| 320 | ! reference wavelength using |
---|
| 321 | ! Qabs instead of Qext |
---|
| 322 | ! (direct comparison with TES) |
---|
[42] | 323 | |
---|
[520] | 324 | |
---|
[411] | 325 | c Test 1d/3d scavenging |
---|
[520] | 326 | real h2otot(ngridmx) |
---|
| 327 | REAL satu(ngridmx,nlayermx) ! satu ratio for output |
---|
| 328 | REAL zqsat(ngridmx,nlayermx) ! saturation |
---|
[42] | 329 | |
---|
| 330 | REAL time_phys |
---|
| 331 | |
---|
[414] | 332 | ! Added for new NLTE scheme |
---|
| 333 | |
---|
| 334 | real co2vmr_gcm(ngridmx,nlayermx) |
---|
| 335 | real n2vmr_gcm(ngridmx,nlayermx) |
---|
| 336 | real ovmr_gcm(ngridmx,nlayermx) |
---|
| 337 | real covmr_gcm(ngridmx,nlayermx) |
---|
| 338 | |
---|
| 339 | |
---|
[267] | 340 | c Variables for PBL |
---|
[544] | 341 | REAL zz1(ngridmx) |
---|
[185] | 342 | REAL lmax_th_out(ngridmx),zmax_th(ngridmx) |
---|
[528] | 343 | REAL, SAVE :: wstar(ngridmx) |
---|
| 344 | REAL, SAVE :: hfmax_th(ngridmx) |
---|
[185] | 345 | REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx) |
---|
| 346 | REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx) |
---|
| 347 | INTEGER lmax_th(ngridmx) |
---|
| 348 | REAL dtke_th(ngridmx,nlayermx+1) |
---|
[268] | 349 | REAL zcdv(ngridmx), zcdh(ngridmx) |
---|
[267] | 350 | REAL Teta_out(ngridmx),u_out(ngridmx) ! Interpolated teta and u at z_out |
---|
[319] | 351 | REAL z_out ! height of interpolation between z0 and z1 [meters] |
---|
[267] | 352 | REAL ustar(ngridmx),tstar(ngridmx) ! friction velocity and friction potential temp |
---|
[528] | 353 | REAL L_mo(ngridmx),wstarpbl(ngridmx),vhf(ngridmx),vvv(ngridmx) |
---|
[267] | 354 | REAL zu2(ngridmx) |
---|
[42] | 355 | c======================================================================= |
---|
| 356 | |
---|
| 357 | c 1. Initialisation: |
---|
| 358 | c ----------------- |
---|
| 359 | |
---|
| 360 | c 1.1 Initialisation only at first call |
---|
| 361 | c --------------------------------------- |
---|
| 362 | IF (firstcall) THEN |
---|
| 363 | |
---|
| 364 | c variables set to 0 |
---|
| 365 | c ~~~~~~~~~~~~~~~~~~ |
---|
[286] | 366 | aerosol(:,:,:)=0 |
---|
| 367 | dtrad(:,:)=0 |
---|
| 368 | fluxrad(:)=0 |
---|
[42] | 369 | |
---|
[528] | 370 | wstar(:)=0. |
---|
[268] | 371 | |
---|
[42] | 372 | c read startfi |
---|
| 373 | c ~~~~~~~~~~~~ |
---|
[226] | 374 | #ifndef MESOSCALE |
---|
| 375 | ! Read netcdf initial physical parameters. |
---|
| 376 | CALL phyetat0 ("startfi.nc",0,0, |
---|
| 377 | & nsoilmx,nq, |
---|
| 378 | & day_ini,time_phys, |
---|
| 379 | & tsurf,tsoil,emis,q2,qsurf,co2ice) |
---|
| 380 | #else |
---|
| 381 | #include "meso_inc/meso_inc_ini.F" |
---|
| 382 | #endif |
---|
[185] | 383 | |
---|
[42] | 384 | if (pday.ne.day_ini) then |
---|
| 385 | write(*,*) "PHYSIQ: ERROR: bad synchronization between ", |
---|
| 386 | & "physics and dynamics" |
---|
| 387 | write(*,*) "dynamics day: ",pday |
---|
| 388 | write(*,*) "physics day: ",day_ini |
---|
| 389 | stop |
---|
| 390 | endif |
---|
| 391 | |
---|
| 392 | write (*,*) 'In physiq day_ini =', day_ini |
---|
| 393 | |
---|
[286] | 394 | c initialize tracers |
---|
| 395 | c ~~~~~~~~~~~~~~~~~~ |
---|
| 396 | tracerdyn=tracer |
---|
| 397 | IF (tracer) THEN |
---|
| 398 | CALL initracer(qsurf,co2ice) |
---|
| 399 | ENDIF ! end tracer |
---|
| 400 | |
---|
[42] | 401 | c Initialize albedo and orbital calculation |
---|
| 402 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 403 | CALL surfini(ngrid,co2ice,qsurf,albedo) |
---|
| 404 | CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) |
---|
| 405 | |
---|
| 406 | c initialize soil |
---|
| 407 | c ~~~~~~~~~~~~~~~ |
---|
| 408 | IF (callsoil) THEN |
---|
| 409 | CALL soil(ngrid,nsoilmx,firstcall,inertiedat, |
---|
| 410 | s ptimestep,tsurf,tsoil,capcal,fluxgrd) |
---|
| 411 | ELSE |
---|
| 412 | PRINT*, |
---|
| 413 | & 'PHYSIQ WARNING! Thermal conduction in the soil turned off' |
---|
| 414 | DO ig=1,ngrid |
---|
| 415 | capcal(ig)=1.e5 |
---|
| 416 | fluxgrd(ig)=0. |
---|
| 417 | ENDDO |
---|
| 418 | ENDIF |
---|
| 419 | icount=1 |
---|
| 420 | |
---|
[226] | 421 | #ifndef MESOSCALE |
---|
| 422 | c Initialize thermospheric parameters |
---|
| 423 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[42] | 424 | |
---|
[226] | 425 | if (callthermos) call param_read |
---|
| 426 | #endif |
---|
[42] | 427 | c Initialize R and Cp as constant |
---|
| 428 | |
---|
| 429 | if (.not.callthermos .and. .not.photochem) then |
---|
| 430 | do l=1,nlayermx |
---|
| 431 | do ig=1,ngridmx |
---|
| 432 | rnew(ig,l)=r |
---|
| 433 | cpnew(ig,l)=cpp |
---|
| 434 | mmean(ig,l)=mugaz |
---|
| 435 | enddo |
---|
| 436 | enddo |
---|
| 437 | endif |
---|
| 438 | |
---|
[414] | 439 | if(callnlte.and.nltemodel.eq.2) call NLTE_leedat |
---|
| 440 | if(callnirco2.and.nircorr.eq.1) call NIR_leedat |
---|
| 441 | |
---|
[42] | 442 | IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN |
---|
[283] | 443 | write(*,*)"physiq: water_param Surface water ice albedo:", |
---|
| 444 | . albedo_h2o_ice |
---|
[42] | 445 | ENDIF |
---|
| 446 | |
---|
| 447 | ENDIF ! (end of "if firstcall") |
---|
| 448 | |
---|
[414] | 449 | |
---|
[42] | 450 | c --------------------------------------------------- |
---|
| 451 | c 1.2 Initializations done at every physical timestep: |
---|
| 452 | c --------------------------------------------------- |
---|
| 453 | c |
---|
| 454 | IF (ngrid.NE.ngridmx) THEN |
---|
| 455 | PRINT*,'STOP in PHYSIQ' |
---|
| 456 | PRINT*,'Probleme de dimensions :' |
---|
| 457 | PRINT*,'ngrid = ',ngrid |
---|
| 458 | PRINT*,'ngridmx = ',ngridmx |
---|
| 459 | STOP |
---|
| 460 | ENDIF |
---|
| 461 | |
---|
| 462 | c Initialize various variables |
---|
| 463 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[286] | 464 | pdv(:,:)=0 |
---|
| 465 | pdu(:,:)=0 |
---|
| 466 | pdt(:,:)=0 |
---|
| 467 | pdq(:,:,:)=0 |
---|
| 468 | pdpsrf(:)=0 |
---|
| 469 | zflubid(:)=0 |
---|
| 470 | zdtsurf(:)=0 |
---|
| 471 | dqsurf(:,:)=0 |
---|
[42] | 472 | igout=ngrid/2+1 |
---|
| 473 | |
---|
| 474 | |
---|
| 475 | zday=pday+ptime ! compute time, in sols (and fraction thereof) |
---|
| 476 | |
---|
| 477 | c Compute Solar Longitude (Ls) : |
---|
| 478 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 479 | if (season) then |
---|
| 480 | call solarlong(zday,zls) |
---|
| 481 | else |
---|
| 482 | call solarlong(float(day_ini),zls) |
---|
| 483 | end if |
---|
| 484 | |
---|
| 485 | c Compute geopotential at interlayers |
---|
| 486 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 487 | c ponderation des altitudes au niveau des couches en dp/p |
---|
| 488 | |
---|
| 489 | DO l=1,nlayer |
---|
| 490 | DO ig=1,ngrid |
---|
| 491 | zzlay(ig,l)=pphi(ig,l)/g |
---|
| 492 | ENDDO |
---|
| 493 | ENDDO |
---|
| 494 | DO ig=1,ngrid |
---|
| 495 | zzlev(ig,1)=0. |
---|
| 496 | zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km... |
---|
| 497 | ENDDO |
---|
| 498 | DO l=2,nlayer |
---|
| 499 | DO ig=1,ngrid |
---|
| 500 | z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) |
---|
| 501 | z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) |
---|
| 502 | zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) |
---|
| 503 | ENDDO |
---|
| 504 | ENDDO |
---|
| 505 | |
---|
| 506 | |
---|
| 507 | ! Potential temperature calculation not the same in physiq and dynamic |
---|
| 508 | |
---|
| 509 | c Compute potential temperature |
---|
| 510 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 511 | DO l=1,nlayer |
---|
| 512 | DO ig=1,ngrid |
---|
| 513 | zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp |
---|
| 514 | zh(ig,l)=pt(ig,l)/zpopsk(ig,l) |
---|
| 515 | ENDDO |
---|
| 516 | ENDDO |
---|
| 517 | |
---|
[226] | 518 | #ifndef MESOSCALE |
---|
| 519 | c----------------------------------------------------------------------- |
---|
| 520 | c 1.2.5 Compute mean mass, cp, and R |
---|
| 521 | c -------------------------------- |
---|
[42] | 522 | |
---|
[226] | 523 | if(photochem.or.callthermos) then |
---|
| 524 | call concentrations(pplay,pt,pdt,pq,pdq,ptimestep) |
---|
| 525 | endif |
---|
| 526 | #endif |
---|
[42] | 527 | c----------------------------------------------------------------------- |
---|
| 528 | c 2. Compute radiative tendencies : |
---|
| 529 | c------------------------------------ |
---|
| 530 | |
---|
| 531 | |
---|
| 532 | IF (callrad) THEN |
---|
| 533 | IF( MOD(icount-1,iradia).EQ.0) THEN |
---|
| 534 | |
---|
| 535 | c Local Solar zenith angle |
---|
| 536 | c ~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 537 | CALL orbite(zls,dist_sol,declin) |
---|
| 538 | |
---|
| 539 | IF(diurnal) THEN |
---|
| 540 | ztim1=SIN(declin) |
---|
| 541 | ztim2=COS(declin)*COS(2.*pi*(zday-.5)) |
---|
| 542 | ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) |
---|
| 543 | |
---|
| 544 | CALL solang(ngrid,sinlon,coslon,sinlat,coslat, |
---|
| 545 | s ztim1,ztim2,ztim3, mu0,fract) |
---|
| 546 | |
---|
| 547 | ELSE |
---|
| 548 | CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad) |
---|
| 549 | ENDIF |
---|
| 550 | |
---|
| 551 | c NLTE cooling from CO2 emission |
---|
| 552 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[414] | 553 | IF(callnlte) then |
---|
| 554 | if(nltemodel.eq.0.or.nltemodel.eq.1) then |
---|
| 555 | CALL nltecool(ngrid,nlayer,nq,pplay,pt,pq,zdtnlte) |
---|
| 556 | else if(nltemodel.eq.2) then |
---|
| 557 | do ig=1,ngrid |
---|
| 558 | do l=1,nlayer |
---|
| 559 | co2vmr_gcm(ig,l)=pq(ig,l,igcm_co2)* |
---|
| 560 | $ mmean(ig,l)/mmol(igcm_co2) |
---|
| 561 | n2vmr_gcm(ig,l)=pq(ig,l,igcm_n2)* |
---|
| 562 | $ mmean(ig,l)/mmol(igcm_n2) |
---|
| 563 | covmr_gcm(ig,l)=pq(ig,l,igcm_co)* |
---|
| 564 | $ mmean(ig,l)/mmol(igcm_co) |
---|
| 565 | ovmr_gcm(ig,l)=pq(ig,l,igcm_o)* |
---|
| 566 | $ mmean(ig,l)/mmol(igcm_o) |
---|
| 567 | enddo |
---|
| 568 | enddo |
---|
| 569 | |
---|
| 570 | CALL NLTEdlvr09_TCOOL(ngrid,nlayer,pplay*9.869e-6, |
---|
| 571 | $ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, |
---|
| 572 | $ ovmr_gcm, zdtnlte ) |
---|
[42] | 573 | |
---|
[414] | 574 | do ig=1,ngrid |
---|
| 575 | do l=1,nlayer |
---|
| 576 | zdtnlte(ig,l)=zdtnlte(ig,l)/86400. |
---|
| 577 | enddo |
---|
| 578 | enddo |
---|
| 579 | endif |
---|
[528] | 580 | else |
---|
| 581 | zdtnlte(:,:)=0. |
---|
[414] | 582 | endif |
---|
[42] | 583 | |
---|
| 584 | c Find number of layers for LTE radiation calculations |
---|
| 585 | IF(MOD(iphysiq*(icount-1),day_step).EQ.0) |
---|
| 586 | & CALL nlthermeq(ngrid,nlayer,pplev,pplay) |
---|
| 587 | |
---|
| 588 | c Note: Dustopacity.F has been transferred to callradite.F |
---|
| 589 | |
---|
| 590 | c Call main radiative transfer scheme |
---|
| 591 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 592 | c Transfer through CO2 (except NIR CO2 absorption) |
---|
| 593 | c and aerosols (dust and water ice) |
---|
| 594 | |
---|
| 595 | c Radiative transfer |
---|
| 596 | c ------------------ |
---|
| 597 | CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, |
---|
| 598 | $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, |
---|
| 599 | $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw, |
---|
[520] | 600 | $ tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice, |
---|
| 601 | $ nuice,co2ice) |
---|
[42] | 602 | |
---|
[234] | 603 | c Outputs for basic check (middle of domain) |
---|
| 604 | c ------------------------------------------ |
---|
[286] | 605 | print*, 'Ls =',zls*180./pi, |
---|
| 606 | & 'check lat lon', lati(igout)*180/pi, |
---|
| 607 | & long(igout)*180/pi |
---|
| 608 | print*, 'tauref(700 Pa) =',tauref(igout), |
---|
| 609 | & ' tau(700 Pa) =',tau(igout,1)*700./pplev(igout,1) |
---|
[42] | 610 | |
---|
[234] | 611 | c --------------------------------------------------------- |
---|
| 612 | c Call slope parameterization for direct and scattered flux |
---|
| 613 | c --------------------------------------------------------- |
---|
| 614 | IF(callslope) THEN |
---|
| 615 | print *, 'Slope scheme is on and computing...' |
---|
| 616 | DO ig=1,ngrid |
---|
| 617 | sl_the = theta_sl(ig) |
---|
| 618 | IF (sl_the .ne. 0.) THEN |
---|
| 619 | ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2) |
---|
| 620 | DO l=1,2 |
---|
| 621 | sl_lct = ptime*24. + 180.*long(ig)/pi/15. |
---|
| 622 | sl_ra = pi*(1.0-sl_lct/12.) |
---|
| 623 | sl_lat = 180.*lati(ig)/pi |
---|
| 624 | sl_tau = tau(ig,1) |
---|
| 625 | sl_alb = albedo(ig,l) |
---|
| 626 | sl_psi = psi_sl(ig) |
---|
| 627 | sl_fl0 = fluxsurf_sw(ig,l) |
---|
| 628 | sl_di0 = 0. |
---|
| 629 | if (mu0(ig) .gt. 0.) then |
---|
| 630 | sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) |
---|
| 631 | sl_di0 = sl_di0*1370./dist_sol/dist_sol |
---|
| 632 | sl_di0 = sl_di0/ztim1 |
---|
| 633 | sl_di0 = fluxsurf_sw(ig,l)*sl_di0 |
---|
| 634 | endif |
---|
| 635 | ! you never know (roundup concern...) |
---|
| 636 | if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 |
---|
| 637 | !!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 638 | CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, |
---|
| 639 | & sl_tau, sl_alb, sl_the, sl_psi, |
---|
| 640 | & sl_di0, sl_fl0, sl_flu ) |
---|
| 641 | !!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 642 | fluxsurf_sw(ig,l) = sl_flu |
---|
| 643 | ENDDO |
---|
| 644 | !!! compute correction on IR flux as well |
---|
| 645 | sky= (1.+cos(pi*theta_sl(ig)/180.))/2. |
---|
| 646 | fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky |
---|
| 647 | ENDIF |
---|
| 648 | ENDDO |
---|
| 649 | ENDIF |
---|
| 650 | |
---|
[42] | 651 | c CO2 near infrared absorption |
---|
| 652 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[286] | 653 | zdtnirco2(:,:)=0 |
---|
[42] | 654 | if (callnirco2) then |
---|
[414] | 655 | call nirco2abs (ngrid,nlayer,pplay,dist_sol,nq,pq, |
---|
[42] | 656 | . mu0,fract,declin, zdtnirco2) |
---|
| 657 | endif |
---|
| 658 | |
---|
| 659 | c Radiative flux from the sky absorbed by the surface (W.m-2) |
---|
| 660 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 661 | DO ig=1,ngrid |
---|
| 662 | fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig) |
---|
| 663 | $ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1)) |
---|
| 664 | $ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2)) |
---|
| 665 | ENDDO |
---|
| 666 | |
---|
| 667 | |
---|
| 668 | c Net atmospheric radiative heating rate (K.s-1) |
---|
| 669 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 670 | IF(callnlte) THEN |
---|
| 671 | CALL blendrad(ngrid, nlayer, pplay, |
---|
| 672 | & zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad) |
---|
| 673 | ELSE |
---|
| 674 | DO l=1,nlayer |
---|
| 675 | DO ig=1,ngrid |
---|
| 676 | dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) |
---|
| 677 | & +zdtnirco2(ig,l) |
---|
| 678 | ENDDO |
---|
| 679 | ENDDO |
---|
| 680 | ENDIF |
---|
| 681 | |
---|
| 682 | ENDIF ! of if(mod(icount-1,iradia).eq.0) |
---|
| 683 | |
---|
| 684 | c Transformation of the radiative tendencies: |
---|
| 685 | c ------------------------------------------- |
---|
| 686 | |
---|
| 687 | c Net radiative surface flux (W.m-2) |
---|
| 688 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 689 | c |
---|
| 690 | DO ig=1,ngrid |
---|
| 691 | zplanck(ig)=tsurf(ig)*tsurf(ig) |
---|
| 692 | zplanck(ig)=emis(ig)* |
---|
| 693 | $ stephan*zplanck(ig)*zplanck(ig) |
---|
| 694 | fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig) |
---|
[234] | 695 | IF(callslope) THEN |
---|
| 696 | sky= (1.+cos(pi*theta_sl(ig)/180.))/2. |
---|
| 697 | fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig) |
---|
| 698 | ENDIF |
---|
[42] | 699 | ENDDO |
---|
| 700 | |
---|
| 701 | DO l=1,nlayer |
---|
| 702 | DO ig=1,ngrid |
---|
| 703 | pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) |
---|
| 704 | ENDDO |
---|
| 705 | ENDDO |
---|
| 706 | |
---|
| 707 | ENDIF ! of IF (callrad) |
---|
| 708 | |
---|
[226] | 709 | c----------------------------------------------------------------------- |
---|
| 710 | c 3. Gravity wave and subgrid scale topography drag : |
---|
| 711 | c ------------------------------------------------- |
---|
[42] | 712 | |
---|
[226] | 713 | |
---|
| 714 | IF(calllott)THEN |
---|
| 715 | |
---|
| 716 | CALL calldrag_noro(ngrid,nlayer,ptimestep, |
---|
| 717 | & pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw) |
---|
| 718 | |
---|
| 719 | DO l=1,nlayer |
---|
| 720 | DO ig=1,ngrid |
---|
| 721 | pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l) |
---|
| 722 | pdu(ig,l)=pdu(ig,l)+zdugw(ig,l) |
---|
| 723 | pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l) |
---|
| 724 | ENDDO |
---|
| 725 | ENDDO |
---|
| 726 | ENDIF |
---|
[234] | 727 | |
---|
[42] | 728 | c----------------------------------------------------------------------- |
---|
| 729 | c 4. Vertical diffusion (turbulent mixing): |
---|
| 730 | c ----------------------------------------- |
---|
[226] | 731 | |
---|
[42] | 732 | IF (calldifv) THEN |
---|
| 733 | |
---|
| 734 | DO ig=1,ngrid |
---|
| 735 | zflubid(ig)=fluxrad(ig)+fluxgrd(ig) |
---|
| 736 | ENDDO |
---|
| 737 | |
---|
[286] | 738 | zdum1(:,:)=0 |
---|
| 739 | zdum2(:,:)=0 |
---|
[42] | 740 | do l=1,nlayer |
---|
| 741 | do ig=1,ngrid |
---|
| 742 | zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) |
---|
| 743 | enddo |
---|
| 744 | enddo |
---|
[226] | 745 | |
---|
[288] | 746 | |
---|
| 747 | #ifdef MESOSCALE |
---|
[289] | 748 | IF (.not.flag_LES) THEN |
---|
[288] | 749 | #endif |
---|
| 750 | c ---------------------- |
---|
[284] | 751 | c Treatment of a special case : using new surface layer (Richardson based) |
---|
| 752 | c without using the thermals in gcm and mesoscale can yield problems in |
---|
| 753 | c weakly unstable situations when winds are near to 0. For those cases, we add |
---|
| 754 | c a unit subgrid gustiness. Remember that thermals should be used we using the |
---|
| 755 | c Richardson based surface layer model. |
---|
[288] | 756 | IF ( .not.calltherm .and. callrichsl ) THEN |
---|
[284] | 757 | DO ig=1, ngridmx |
---|
| 758 | IF (zh(ig,1) .lt. tsurf(ig)) THEN |
---|
[528] | 759 | wstar(ig)=1. |
---|
| 760 | hfmax_th(ig)=0.2 |
---|
| 761 | ELSE |
---|
| 762 | wstar(ig)=0. |
---|
| 763 | hfmax_th(ig)=0. |
---|
| 764 | ENDIF |
---|
[284] | 765 | ENDDO |
---|
| 766 | ENDIF |
---|
[288] | 767 | c ---------------------- |
---|
| 768 | #ifdef MESOSCALE |
---|
[284] | 769 | ENDIF |
---|
| 770 | #endif |
---|
| 771 | |
---|
[544] | 772 | IF (tke_heat_flux .ne. 0.) THEN |
---|
| 773 | zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)* |
---|
| 774 | & (-alog(pplay(:,1)/pplev(:,1))) |
---|
| 775 | pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1) |
---|
| 776 | ENDIF |
---|
[288] | 777 | |
---|
[42] | 778 | c Calling vdif (Martian version WITH CO2 condensation) |
---|
| 779 | CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk, |
---|
| 780 | $ ptimestep,capcal,lwrite, |
---|
| 781 | $ pplay,pplev,zzlay,zzlev,z0, |
---|
| 782 | $ pu,pv,zh,pq,tsurf,emis,qsurf, |
---|
| 783 | $ zdum1,zdum2,zdh,pdq,zflubid, |
---|
| 784 | $ zdudif,zdvdif,zdhdif,zdtsdif,q2, |
---|
[529] | 785 | & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th |
---|
| 786 | #ifdef MESOSCALE |
---|
| 787 | & ,flag_LES |
---|
| 788 | #endif |
---|
| 789 | & ) |
---|
[42] | 790 | |
---|
[529] | 791 | |
---|
[226] | 792 | #ifdef MESOSCALE |
---|
| 793 | #include "meso_inc/meso_inc_les.F" |
---|
| 794 | #endif |
---|
[42] | 795 | DO l=1,nlayer |
---|
| 796 | DO ig=1,ngrid |
---|
| 797 | pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l) |
---|
| 798 | pdu(ig,l)=pdu(ig,l)+zdudif(ig,l) |
---|
| 799 | pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l) |
---|
| 800 | |
---|
| 801 | zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only |
---|
| 802 | |
---|
| 803 | ENDDO |
---|
| 804 | ENDDO |
---|
| 805 | |
---|
[226] | 806 | DO ig=1,ngrid |
---|
| 807 | zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig) |
---|
| 808 | ENDDO |
---|
[42] | 809 | |
---|
| 810 | if (tracer) then |
---|
| 811 | DO iq=1, nq |
---|
| 812 | DO l=1,nlayer |
---|
| 813 | DO ig=1,ngrid |
---|
| 814 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) |
---|
| 815 | ENDDO |
---|
| 816 | ENDDO |
---|
| 817 | ENDDO |
---|
| 818 | DO iq=1, nq |
---|
| 819 | DO ig=1,ngrid |
---|
| 820 | dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq) |
---|
| 821 | ENDDO |
---|
| 822 | ENDDO |
---|
| 823 | end if ! of if (tracer) |
---|
| 824 | |
---|
| 825 | ELSE |
---|
| 826 | DO ig=1,ngrid |
---|
| 827 | zdtsurf(ig)=zdtsurf(ig)+ |
---|
| 828 | s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) |
---|
| 829 | ENDDO |
---|
[226] | 830 | #ifdef MESOSCALE |
---|
[42] | 831 | IF (flag_LES) THEN |
---|
| 832 | write(*,*) 'LES mode !' |
---|
| 833 | write(*,*) 'Please set calldifv to T in callphys.def' |
---|
| 834 | STOP |
---|
| 835 | ENDIF |
---|
[226] | 836 | #endif |
---|
[42] | 837 | ENDIF ! of IF (calldifv) |
---|
| 838 | |
---|
[162] | 839 | c----------------------------------------------------------------------- |
---|
| 840 | c TEST. Thermals : |
---|
| 841 | c HIGHLY EXPERIMENTAL, BEWARE !! |
---|
| 842 | c ----------------------------- |
---|
| 843 | |
---|
| 844 | if(calltherm) then |
---|
| 845 | |
---|
[185] | 846 | call calltherm_interface(firstcall, |
---|
[162] | 847 | $ long,lati,zzlev,zzlay, |
---|
| 848 | $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, |
---|
[185] | 849 | $ pplay,pplev,pphi,zpopsk, |
---|
| 850 | $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th, |
---|
[528] | 851 | $ dtke_th,hfmax_th,wstar) |
---|
[226] | 852 | |
---|
[162] | 853 | DO l=1,nlayer |
---|
| 854 | DO ig=1,ngrid |
---|
| 855 | pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l) |
---|
| 856 | pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l) |
---|
| 857 | pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l) |
---|
| 858 | q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep |
---|
| 859 | ENDDO |
---|
| 860 | ENDDO |
---|
| 861 | |
---|
| 862 | DO ig=1,ngrid |
---|
| 863 | q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep |
---|
| 864 | ENDDO |
---|
| 865 | |
---|
| 866 | if (tracer) then |
---|
| 867 | DO iq=1,nq |
---|
| 868 | DO l=1,nlayer |
---|
| 869 | DO ig=1,ngrid |
---|
| 870 | pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq) |
---|
| 871 | ENDDO |
---|
| 872 | ENDDO |
---|
| 873 | ENDDO |
---|
| 874 | endif |
---|
[42] | 875 | |
---|
[277] | 876 | lmax_th_out(:)=real(lmax_th(:)) |
---|
| 877 | |
---|
[162] | 878 | else !of if calltherm |
---|
| 879 | lmax_th(:)=0 |
---|
[528] | 880 | wstar(:)=0. |
---|
| 881 | hfmax_th(:)=0. |
---|
[277] | 882 | lmax_th_out(:)=0. |
---|
[162] | 883 | end if |
---|
| 884 | |
---|
[42] | 885 | c----------------------------------------------------------------------- |
---|
| 886 | c 5. Dry convective adjustment: |
---|
| 887 | c ----------------------------- |
---|
| 888 | |
---|
| 889 | IF(calladj) THEN |
---|
| 890 | |
---|
| 891 | DO l=1,nlayer |
---|
| 892 | DO ig=1,ngrid |
---|
| 893 | zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) |
---|
| 894 | ENDDO |
---|
| 895 | ENDDO |
---|
[286] | 896 | zduadj(:,:)=0 |
---|
| 897 | zdvadj(:,:)=0 |
---|
| 898 | zdhadj(:,:)=0 |
---|
[42] | 899 | |
---|
| 900 | CALL convadj(ngrid,nlayer,nq,ptimestep, |
---|
[162] | 901 | $ pplay,pplev,zpopsk,lmax_th, |
---|
[42] | 902 | $ pu,pv,zh,pq, |
---|
| 903 | $ pdu,pdv,zdh,pdq, |
---|
| 904 | $ zduadj,zdvadj,zdhadj, |
---|
| 905 | $ zdqadj) |
---|
| 906 | |
---|
[162] | 907 | |
---|
[42] | 908 | DO l=1,nlayer |
---|
| 909 | DO ig=1,ngrid |
---|
| 910 | pdu(ig,l)=pdu(ig,l)+zduadj(ig,l) |
---|
| 911 | pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l) |
---|
| 912 | pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l) |
---|
| 913 | |
---|
| 914 | zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only |
---|
| 915 | ENDDO |
---|
| 916 | ENDDO |
---|
| 917 | |
---|
| 918 | if(tracer) then |
---|
| 919 | DO iq=1, nq |
---|
| 920 | DO l=1,nlayer |
---|
| 921 | DO ig=1,ngrid |
---|
| 922 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) |
---|
| 923 | ENDDO |
---|
| 924 | ENDDO |
---|
| 925 | ENDDO |
---|
| 926 | end if |
---|
| 927 | ENDIF ! of IF(calladj) |
---|
| 928 | |
---|
| 929 | c----------------------------------------------------------------------- |
---|
| 930 | c 6. Carbon dioxide condensation-sublimation: |
---|
| 931 | c ------------------------------------------- |
---|
| 932 | |
---|
[528] | 933 | IF (tituscap) THEN |
---|
| 934 | !!! get the actual co2 seasonal cap from Titus observations |
---|
| 935 | CALL geticecover( ngrid, 180.*zls/pi, |
---|
[333] | 936 | . 180.*long/pi, 180.*lati/pi, co2ice ) |
---|
[528] | 937 | co2ice = co2ice * 10000. |
---|
| 938 | ENDIF |
---|
[332] | 939 | |
---|
[42] | 940 | IF (callcond) THEN |
---|
[86] | 941 | CALL newcondens(ngrid,nlayer,nq,ptimestep, |
---|
[42] | 942 | $ capcal,pplay,pplev,tsurf,pt, |
---|
| 943 | $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, |
---|
| 944 | $ co2ice,albedo,emis, |
---|
| 945 | $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, |
---|
[226] | 946 | $ fluxsurf_sw,zls) |
---|
[42] | 947 | |
---|
| 948 | DO l=1,nlayer |
---|
| 949 | DO ig=1,ngrid |
---|
| 950 | pdt(ig,l)=pdt(ig,l)+zdtc(ig,l) |
---|
| 951 | pdv(ig,l)=pdv(ig,l)+zdvc(ig,l) |
---|
| 952 | pdu(ig,l)=pdu(ig,l)+zduc(ig,l) |
---|
| 953 | ENDDO |
---|
| 954 | ENDDO |
---|
| 955 | DO ig=1,ngrid |
---|
| 956 | zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig) |
---|
| 957 | ENDDO |
---|
| 958 | |
---|
| 959 | IF (tracer) THEN |
---|
| 960 | DO iq=1, nq |
---|
| 961 | DO l=1,nlayer |
---|
| 962 | DO ig=1,ngrid |
---|
| 963 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) |
---|
| 964 | ENDDO |
---|
| 965 | ENDDO |
---|
| 966 | ENDDO |
---|
| 967 | ENDIF ! of IF (tracer) |
---|
| 968 | |
---|
| 969 | ENDIF ! of IF (callcond) |
---|
| 970 | |
---|
| 971 | c----------------------------------------------------------------------- |
---|
| 972 | c 7. Specific parameterizations for tracers |
---|
| 973 | c: ----------------------------------------- |
---|
| 974 | |
---|
| 975 | if (tracer) then |
---|
| 976 | |
---|
| 977 | c 7a. Water and ice |
---|
| 978 | c --------------- |
---|
| 979 | |
---|
| 980 | c --------------------------------------- |
---|
| 981 | c Water ice condensation in the atmosphere |
---|
| 982 | c ---------------------------------------- |
---|
| 983 | IF (water) THEN |
---|
| 984 | |
---|
| 985 | call watercloud(ngrid,nlayer,ptimestep, |
---|
| 986 | & pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, |
---|
| 987 | & pq,pdq,zdqcloud,zdqscloud,zdtcloud, |
---|
[358] | 988 | & nq,tau,tauscaling,rdust,rice,nuice, |
---|
| 989 | & rsedcloud,rhocloud) |
---|
[42] | 990 | if (activice) then |
---|
| 991 | c Temperature variation due to latent heat release |
---|
| 992 | DO l=1,nlayer |
---|
| 993 | DO ig=1,ngrid |
---|
| 994 | pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l) |
---|
| 995 | ENDDO |
---|
| 996 | ENDDO |
---|
| 997 | endif |
---|
[528] | 998 | |
---|
[42] | 999 | ! increment water vapour and ice atmospheric tracers tendencies |
---|
| 1000 | IF (water) THEN |
---|
[358] | 1001 | DO l=1,nlayer |
---|
[42] | 1002 | DO ig=1,ngrid |
---|
[358] | 1003 | pdq(ig,l,igcm_h2o_vap)= |
---|
| 1004 | & pdq(ig,l,igcm_h2o_vap)+ |
---|
| 1005 | & zdqcloud(ig,l,igcm_h2o_vap) |
---|
| 1006 | pdq(ig,l,igcm_h2o_ice)= |
---|
| 1007 | & pdq(ig,l,igcm_h2o_ice)+ |
---|
| 1008 | & zdqcloud(ig,l,igcm_h2o_ice) |
---|
[520] | 1009 | if (((pq(ig,l,igcm_h2o_ice) + |
---|
| 1010 | & pdq(ig,l,igcm_h2o_ice)*ptimestep)) .le. 0) |
---|
| 1011 | & then |
---|
| 1012 | pdq(ig,l,igcm_h2o_ice) = |
---|
| 1013 | & - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30 |
---|
| 1014 | endif |
---|
[358] | 1015 | IF (scavenging) THEN |
---|
| 1016 | pdq(ig,l,igcm_ccn_mass)= |
---|
| 1017 | & pdq(ig,l,igcm_ccn_mass)+ |
---|
| 1018 | & zdqcloud(ig,l,igcm_ccn_mass) |
---|
| 1019 | pdq(ig,l,igcm_ccn_number)= |
---|
| 1020 | & pdq(ig,l,igcm_ccn_number)+ |
---|
| 1021 | & zdqcloud(ig,l,igcm_ccn_number) |
---|
[520] | 1022 | pdq(ig,l,igcm_dust_mass)= |
---|
| 1023 | & pdq(ig,l,igcm_dust_mass)+ |
---|
| 1024 | & zdqcloud(ig,l,igcm_dust_mass) |
---|
| 1025 | pdq(ig,l,igcm_dust_number)= |
---|
| 1026 | & pdq(ig,l,igcm_dust_number)+ |
---|
| 1027 | & zdqcloud(ig,l,igcm_dust_number) |
---|
[411] | 1028 | !!!!!!!!!!!!!!!!!!!!! We need to check that we have Nccn & Ndust > 0 |
---|
| 1029 | !!!!!!!!!!!!!!!!!!!!! This is due to single precision rounding problems |
---|
[358] | 1030 | if (((pq(ig,l,igcm_ccn_number) + |
---|
| 1031 | & pdq(ig,l,igcm_ccn_number)*ptimestep)) .le. 0) |
---|
| 1032 | & then |
---|
| 1033 | pdq(ig,l,igcm_ccn_number) = |
---|
[411] | 1034 | & - pq(ig,l,igcm_ccn_number)/ptimestep + 1e-30 |
---|
[520] | 1035 | pdq(ig,l,igcm_ccn_mass) = |
---|
| 1036 | & - pq(ig,l,igcm_ccn_mass)/ptimestep + 1e-30 |
---|
[358] | 1037 | endif |
---|
[411] | 1038 | if (((pq(ig,l,igcm_dust_number) + |
---|
| 1039 | & pdq(ig,l,igcm_dust_number)*ptimestep)) .le. 0) |
---|
| 1040 | & then |
---|
| 1041 | pdq(ig,l,igcm_dust_number) = |
---|
| 1042 | & - pq(ig,l,igcm_dust_number)/ptimestep + 1e-30 |
---|
[520] | 1043 | pdq(ig,l,igcm_dust_mass) = |
---|
| 1044 | & - pq(ig,l,igcm_dust_mass)/ptimestep + 1e-30 |
---|
[411] | 1045 | endif |
---|
[358] | 1046 | !!!!!!!!!!!!!!!!!!!!! |
---|
| 1047 | !!!!!!!!!!!!!!!!!!!!! |
---|
| 1048 | ENDIF |
---|
[42] | 1049 | ENDDO |
---|
| 1050 | ENDDO |
---|
| 1051 | ENDIF ! of IF (water) THEN |
---|
[411] | 1052 | |
---|
[42] | 1053 | ! Increment water ice surface tracer tendency |
---|
| 1054 | DO ig=1,ngrid |
---|
| 1055 | dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+ |
---|
| 1056 | & zdqscloud(ig,igcm_h2o_ice) |
---|
| 1057 | ENDDO |
---|
| 1058 | |
---|
| 1059 | END IF ! of IF (water) |
---|
| 1060 | |
---|
[520] | 1061 | |
---|
| 1062 | c 7b. Chemical species |
---|
| 1063 | c ------------------ |
---|
| 1064 | |
---|
| 1065 | #ifndef MESOSCALE |
---|
| 1066 | c -------------- |
---|
| 1067 | c photochemistry : |
---|
| 1068 | c -------------- |
---|
| 1069 | IF (photochem .or. thermochem) then |
---|
| 1070 | !NB: Photochemistry includes condensation of H2O2 |
---|
| 1071 | PRINT*, 'SURFDUST,SURFICE TO BE IMPLEMENTED. YAAAAAARG.' |
---|
| 1072 | call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, |
---|
| 1073 | $ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, |
---|
| 1074 | $ zdqcloud,zdqscloud,tauref,co2ice, |
---|
| 1075 | $ pu,pdu,pv,pdv,surfdust,surfice) |
---|
| 1076 | |
---|
| 1077 | ! increment values of tracers: |
---|
| 1078 | DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry |
---|
| 1079 | ! tracers is zero anyways |
---|
| 1080 | DO l=1,nlayer |
---|
| 1081 | DO ig=1,ngrid |
---|
| 1082 | pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) |
---|
| 1083 | ENDDO |
---|
| 1084 | ENDDO |
---|
| 1085 | ENDDO ! of DO iq=1,nq |
---|
| 1086 | ! add condensation tendency for H2O2 |
---|
| 1087 | if (igcm_h2o2.ne.0) then |
---|
| 1088 | DO l=1,nlayer |
---|
| 1089 | DO ig=1,ngrid |
---|
| 1090 | pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2) |
---|
| 1091 | & +zdqcloud(ig,l,igcm_h2o2) |
---|
| 1092 | ENDDO |
---|
| 1093 | ENDDO |
---|
| 1094 | endif |
---|
| 1095 | |
---|
| 1096 | ! increment surface values of tracers: |
---|
| 1097 | DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry |
---|
| 1098 | ! tracers is zero anyways |
---|
| 1099 | DO ig=1,ngrid |
---|
| 1100 | dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) |
---|
| 1101 | ENDDO |
---|
| 1102 | ENDDO ! of DO iq=1,nq |
---|
| 1103 | ! add condensation tendency for H2O2 |
---|
| 1104 | if (igcm_h2o2.ne.0) then |
---|
| 1105 | DO ig=1,ngrid |
---|
| 1106 | dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2) |
---|
| 1107 | & +zdqscloud(ig,igcm_h2o2) |
---|
| 1108 | ENDDO |
---|
| 1109 | endif |
---|
| 1110 | |
---|
| 1111 | END IF ! of IF (photochem.or.thermochem) |
---|
| 1112 | #endif |
---|
| 1113 | |
---|
| 1114 | c 7c. Aerosol particles |
---|
[42] | 1115 | c ------------------- |
---|
| 1116 | |
---|
| 1117 | c ---------- |
---|
| 1118 | c Dust devil : |
---|
| 1119 | c ---------- |
---|
| 1120 | IF(callddevil) then |
---|
| 1121 | call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2, |
---|
| 1122 | & zdqdev,zdqsdev) |
---|
| 1123 | |
---|
| 1124 | if (dustbin.ge.1) then |
---|
| 1125 | do iq=1,nq |
---|
| 1126 | DO l=1,nlayer |
---|
| 1127 | DO ig=1,ngrid |
---|
| 1128 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq) |
---|
| 1129 | ENDDO |
---|
| 1130 | ENDDO |
---|
| 1131 | enddo |
---|
| 1132 | do iq=1,nq |
---|
| 1133 | DO ig=1,ngrid |
---|
| 1134 | dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq) |
---|
| 1135 | ENDDO |
---|
| 1136 | enddo |
---|
| 1137 | endif ! of if (dustbin.ge.1) |
---|
| 1138 | |
---|
| 1139 | END IF ! of IF (callddevil) |
---|
| 1140 | |
---|
| 1141 | c ------------- |
---|
| 1142 | c Sedimentation : acts also on water ice |
---|
| 1143 | c ------------- |
---|
| 1144 | IF (sedimentation) THEN |
---|
| 1145 | !call zerophys(ngrid*nlayer*nq, zdqsed) |
---|
| 1146 | zdqsed(1:ngrid,1:nlayer,1:nq)=0 |
---|
| 1147 | !call zerophys(ngrid*nq, zdqssed) |
---|
| 1148 | zdqssed(1:ngrid,1:nq)=0 |
---|
| 1149 | |
---|
| 1150 | call callsedim(ngrid,nlayer, ptimestep, |
---|
[411] | 1151 | & pplev,zzlev, zzlay, pt, rdust, rice, |
---|
[358] | 1152 | & rsedcloud,rhocloud, |
---|
[411] | 1153 | & pq, pdq, zdqsed, zdqssed,nq, |
---|
| 1154 | & tau,tauscaling) |
---|
[520] | 1155 | |
---|
| 1156 | |
---|
| 1157 | !print*, 'h2o_ice zdqsed ds physiq', zdqsed(1,:,igcm_h2o_ice) |
---|
[42] | 1158 | |
---|
| 1159 | DO iq=1, nq |
---|
| 1160 | DO l=1,nlayer |
---|
| 1161 | DO ig=1,ngrid |
---|
| 1162 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq) |
---|
| 1163 | ENDDO |
---|
| 1164 | ENDDO |
---|
| 1165 | ENDDO |
---|
| 1166 | DO iq=1, nq |
---|
| 1167 | DO ig=1,ngrid |
---|
| 1168 | dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq) |
---|
| 1169 | ENDDO |
---|
| 1170 | ENDDO |
---|
| 1171 | END IF ! of IF (sedimentation) |
---|
[520] | 1172 | |
---|
[42] | 1173 | |
---|
[411] | 1174 | |
---|
[42] | 1175 | c 7d. Updates |
---|
| 1176 | c --------- |
---|
| 1177 | |
---|
| 1178 | DO iq=1, nq |
---|
| 1179 | DO ig=1,ngrid |
---|
| 1180 | |
---|
| 1181 | c --------------------------------- |
---|
| 1182 | c Updating tracer budget on surface |
---|
| 1183 | c --------------------------------- |
---|
| 1184 | qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq) |
---|
| 1185 | |
---|
| 1186 | ENDDO ! (ig) |
---|
| 1187 | ENDDO ! (iq) |
---|
| 1188 | |
---|
| 1189 | endif ! of if (tracer) |
---|
| 1190 | |
---|
[226] | 1191 | #ifndef MESOSCALE |
---|
| 1192 | c----------------------------------------------------------------------- |
---|
| 1193 | c 8. THERMOSPHERE CALCULATION |
---|
| 1194 | c----------------------------------------------------------------------- |
---|
[42] | 1195 | |
---|
[226] | 1196 | if (callthermos) then |
---|
| 1197 | call thermosphere(pplev,pplay,dist_sol, |
---|
| 1198 | $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay, |
---|
| 1199 | & pt,pq,pu,pv,pdt,pdq, |
---|
| 1200 | $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff) |
---|
| 1201 | |
---|
| 1202 | DO l=1,nlayer |
---|
| 1203 | DO ig=1,ngrid |
---|
| 1204 | dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l) |
---|
| 1205 | pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l) |
---|
| 1206 | & +zdteuv(ig,l) |
---|
| 1207 | pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l) |
---|
| 1208 | pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l) |
---|
| 1209 | DO iq=1, nq |
---|
| 1210 | pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq) |
---|
| 1211 | ENDDO |
---|
| 1212 | ENDDO |
---|
| 1213 | ENDDO |
---|
| 1214 | |
---|
| 1215 | endif ! of if (callthermos) |
---|
| 1216 | #endif |
---|
| 1217 | |
---|
[42] | 1218 | c----------------------------------------------------------------------- |
---|
| 1219 | c 9. Surface and sub-surface soil temperature |
---|
| 1220 | c----------------------------------------------------------------------- |
---|
| 1221 | c |
---|
| 1222 | c |
---|
| 1223 | c 9.1 Increment Surface temperature: |
---|
| 1224 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1225 | |
---|
| 1226 | DO ig=1,ngrid |
---|
| 1227 | tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) |
---|
| 1228 | ENDDO |
---|
| 1229 | |
---|
| 1230 | c Prescribe a cold trap at south pole (except at high obliquity !!) |
---|
| 1231 | c Temperature at the surface is set there to be the temperature |
---|
| 1232 | c corresponding to equilibrium temperature between phases of CO2 |
---|
| 1233 | |
---|
[411] | 1234 | |
---|
[42] | 1235 | IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN |
---|
[226] | 1236 | #ifndef MESOSCALE |
---|
| 1237 | if (caps.and.(obliquit.lt.27.)) then |
---|
| 1238 | ! NB: Updated surface pressure, at grid point 'ngrid', is |
---|
| 1239 | ! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep |
---|
| 1240 | tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095* |
---|
| 1241 | & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) |
---|
| 1242 | endif |
---|
| 1243 | #endif |
---|
[42] | 1244 | c ------------------------------------------------------------- |
---|
[283] | 1245 | c Change of surface albedo in case of ground frost |
---|
[42] | 1246 | c everywhere except on the north permanent cap and in regions |
---|
| 1247 | c covered by dry ice. |
---|
| 1248 | c ALWAYS PLACE these lines after newcondens !!! |
---|
| 1249 | c ------------------------------------------------------------- |
---|
| 1250 | do ig=1,ngrid |
---|
| 1251 | if ((co2ice(ig).eq.0).and. |
---|
[283] | 1252 | & (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then |
---|
| 1253 | albedo(ig,1) = albedo_h2o_ice |
---|
| 1254 | albedo(ig,2) = albedo_h2o_ice |
---|
| 1255 | c write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice) |
---|
| 1256 | c write(*,*) "physiq.F frost :" |
---|
| 1257 | c & ,lati(ig)*180./pi, long(ig)*180./pi |
---|
[42] | 1258 | endif |
---|
| 1259 | enddo ! of do ig=1,ngrid |
---|
| 1260 | ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1)) |
---|
| 1261 | |
---|
[528] | 1262 | c |
---|
[42] | 1263 | c 9.2 Compute soil temperatures and subsurface heat flux: |
---|
| 1264 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1265 | IF (callsoil) THEN |
---|
| 1266 | CALL soil(ngrid,nsoilmx,.false.,inertiedat, |
---|
| 1267 | & ptimestep,tsurf,tsoil,capcal,fluxgrd) |
---|
| 1268 | ENDIF |
---|
[411] | 1269 | |
---|
[42] | 1270 | |
---|
| 1271 | c----------------------------------------------------------------------- |
---|
| 1272 | c 10. Write output files |
---|
| 1273 | c ---------------------- |
---|
| 1274 | |
---|
| 1275 | c ------------------------------- |
---|
| 1276 | c Dynamical fields incrementation |
---|
| 1277 | c ------------------------------- |
---|
| 1278 | c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics) |
---|
| 1279 | ! temperature, zonal and meridional wind |
---|
| 1280 | DO l=1,nlayer |
---|
| 1281 | DO ig=1,ngrid |
---|
| 1282 | zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep |
---|
| 1283 | zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep |
---|
| 1284 | zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep |
---|
| 1285 | ENDDO |
---|
| 1286 | ENDDO |
---|
| 1287 | |
---|
| 1288 | ! tracers |
---|
| 1289 | DO iq=1, nq |
---|
| 1290 | DO l=1,nlayer |
---|
| 1291 | DO ig=1,ngrid |
---|
| 1292 | zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep |
---|
| 1293 | ENDDO |
---|
| 1294 | ENDDO |
---|
| 1295 | ENDDO |
---|
| 1296 | |
---|
| 1297 | ! surface pressure |
---|
| 1298 | DO ig=1,ngrid |
---|
| 1299 | ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep |
---|
| 1300 | ENDDO |
---|
| 1301 | |
---|
| 1302 | ! pressure |
---|
| 1303 | DO l=1,nlayer |
---|
| 1304 | DO ig=1,ngrid |
---|
| 1305 | zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig) |
---|
| 1306 | zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig) |
---|
| 1307 | ENDDO |
---|
| 1308 | ENDDO |
---|
| 1309 | |
---|
| 1310 | ! Density |
---|
| 1311 | DO l=1,nlayer |
---|
| 1312 | DO ig=1,ngrid |
---|
| 1313 | rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l)) |
---|
| 1314 | ENDDO |
---|
| 1315 | ENDDO |
---|
| 1316 | |
---|
[269] | 1317 | ! Potential Temperature |
---|
| 1318 | |
---|
| 1319 | DO ig=1,ngridmx |
---|
| 1320 | DO l=1,nlayermx |
---|
[528] | 1321 | zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp |
---|
[269] | 1322 | ENDDO |
---|
| 1323 | ENDDO |
---|
| 1324 | |
---|
| 1325 | |
---|
[226] | 1326 | c Compute surface stress : (NB: z0 is a common in surfdat.h) |
---|
[42] | 1327 | c DO ig=1,ngrid |
---|
[226] | 1328 | c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2 |
---|
[42] | 1329 | c zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2) |
---|
| 1330 | c ENDDO |
---|
| 1331 | |
---|
| 1332 | c Sum of fluxes in solar spectral bands (for output only) |
---|
| 1333 | DO ig=1,ngrid |
---|
[226] | 1334 | fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2) |
---|
| 1335 | fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2) |
---|
[42] | 1336 | ENDDO |
---|
| 1337 | c ******* TEST ****************************************************** |
---|
| 1338 | ztim1 = 999 |
---|
| 1339 | DO l=1,nlayer |
---|
| 1340 | DO ig=1,ngrid |
---|
| 1341 | if (pt(ig,l).lt.ztim1) then |
---|
| 1342 | ztim1 = pt(ig,l) |
---|
| 1343 | igmin = ig |
---|
| 1344 | lmin = l |
---|
| 1345 | end if |
---|
| 1346 | ENDDO |
---|
| 1347 | ENDDO |
---|
[234] | 1348 | if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then |
---|
[42] | 1349 | write(*,*) 'PHYSIQ: stability WARNING :' |
---|
| 1350 | write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin), |
---|
| 1351 | & 'ig l =', igmin, lmin |
---|
| 1352 | end if |
---|
| 1353 | c ******************************************************************* |
---|
| 1354 | |
---|
| 1355 | c --------------------- |
---|
| 1356 | c Outputs to the screen |
---|
| 1357 | c --------------------- |
---|
| 1358 | |
---|
| 1359 | IF (lwrite) THEN |
---|
| 1360 | PRINT*,'Global diagnostics for the physics' |
---|
| 1361 | PRINT*,'Variables and their increments x and dx/dt * dt' |
---|
| 1362 | WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps' |
---|
| 1363 | WRITE(*,'(2f10.5,2f15.5)') |
---|
| 1364 | s tsurf(igout),zdtsurf(igout)*ptimestep, |
---|
| 1365 | s pplev(igout,1),pdpsrf(igout)*ptimestep |
---|
| 1366 | WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT' |
---|
| 1367 | WRITE(*,'(i4,6f10.5)') (l, |
---|
| 1368 | s pu(igout,l),pdu(igout,l)*ptimestep, |
---|
| 1369 | s pv(igout,l),pdv(igout,l)*ptimestep, |
---|
| 1370 | s pt(igout,l),pdt(igout,l)*ptimestep, |
---|
| 1371 | s l=1,nlayer) |
---|
| 1372 | ENDIF ! of IF (lwrite) |
---|
| 1373 | |
---|
| 1374 | IF (ngrid.NE.1) THEN |
---|
| 1375 | |
---|
[226] | 1376 | #ifndef MESOSCALE |
---|
[42] | 1377 | c ------------------------------------------------------------------- |
---|
| 1378 | c Writing NetCDF file "RESTARTFI" at the end of the run |
---|
| 1379 | c ------------------------------------------------------------------- |
---|
| 1380 | c Note: 'restartfi' is stored just before dynamics are stored |
---|
| 1381 | c in 'restart'. Between now and the writting of 'restart', |
---|
| 1382 | c there will have been the itau=itau+1 instruction and |
---|
| 1383 | c a reset of 'time' (lastacll = .true. when itau+1= itaufin) |
---|
| 1384 | c thus we store for time=time+dtvr |
---|
| 1385 | |
---|
[226] | 1386 | IF(lastcall) THEN |
---|
| 1387 | ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) |
---|
| 1388 | write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin |
---|
| 1389 | call physdem1("restartfi.nc",long,lati,nsoilmx,nq, |
---|
| 1390 | . ptimestep,pday, |
---|
| 1391 | . ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf, |
---|
| 1392 | . area,albedodat,inertiedat,zmea,zstd,zsig, |
---|
| 1393 | . zgam,zthe) |
---|
| 1394 | ENDIF |
---|
| 1395 | #endif |
---|
[42] | 1396 | |
---|
| 1397 | c ------------------------------------------------------------------- |
---|
| 1398 | c Calculation of diagnostic variables written in both stats and |
---|
| 1399 | c diagfi files |
---|
| 1400 | c ------------------------------------------------------------------- |
---|
| 1401 | |
---|
| 1402 | if (tracer) then |
---|
| 1403 | if (water) then |
---|
[520] | 1404 | |
---|
| 1405 | if (scavenging) then |
---|
| 1406 | ccntot(:)= 0 |
---|
| 1407 | do ig=1,ngrid |
---|
| 1408 | do l=1,nlayermx |
---|
| 1409 | ccntot(ig) = ccntot(ig) + |
---|
| 1410 | & zq(ig,l,igcm_ccn_number)*tauscaling(ig) |
---|
| 1411 | & *(pplev(ig,l) - pplev(ig,l+1)) / g |
---|
| 1412 | enddo |
---|
| 1413 | enddo |
---|
| 1414 | endif |
---|
[42] | 1415 | |
---|
[286] | 1416 | mtot(:)=0 |
---|
| 1417 | icetot(:)=0 |
---|
| 1418 | rave(:)=0 |
---|
| 1419 | tauTES(:)=0 |
---|
[42] | 1420 | do ig=1,ngrid |
---|
| 1421 | do l=1,nlayermx |
---|
| 1422 | mtot(ig) = mtot(ig) + |
---|
| 1423 | & zq(ig,l,igcm_h2o_vap) * |
---|
| 1424 | & (pplev(ig,l) - pplev(ig,l+1)) / g |
---|
| 1425 | icetot(ig) = icetot(ig) + |
---|
| 1426 | & zq(ig,l,igcm_h2o_ice) * |
---|
| 1427 | & (pplev(ig,l) - pplev(ig,l+1)) / g |
---|
| 1428 | rave(ig) = rave(ig) + |
---|
| 1429 | & zq(ig,l,igcm_h2o_ice) * |
---|
| 1430 | & (pplev(ig,l) - pplev(ig,l+1)) / g * |
---|
| 1431 | & rice(ig,l) * (1.+nuice_ref) |
---|
| 1432 | c Computing abs optical depth at 825 cm-1 in each |
---|
| 1433 | c layer to simulate NEW TES retrieval |
---|
| 1434 | Qabsice = min( |
---|
| 1435 | & max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2 |
---|
| 1436 | & ) |
---|
| 1437 | opTES(ig,l)= 0.75 * Qabsice * |
---|
| 1438 | & zq(ig,l,igcm_h2o_ice) * |
---|
| 1439 | & (pplev(ig,l) - pplev(ig,l+1)) / g |
---|
| 1440 | & / (rho_ice * rice(ig,l) * (1.+nuice_ref)) |
---|
| 1441 | tauTES(ig)=tauTES(ig)+ opTES(ig,l) |
---|
| 1442 | enddo |
---|
| 1443 | rave(ig)=rave(ig)/max(icetot(ig),1.e-30) |
---|
| 1444 | if (icetot(ig)*1e3.lt.0.01) rave(ig)=0. |
---|
| 1445 | enddo |
---|
| 1446 | |
---|
| 1447 | endif ! of if (water) |
---|
| 1448 | endif ! of if (tracer) |
---|
| 1449 | |
---|
| 1450 | c ----------------------------------------------------------------- |
---|
| 1451 | c WSTATS: Saving statistics |
---|
| 1452 | c ----------------------------------------------------------------- |
---|
| 1453 | c ("stats" stores and accumulates 8 key variables in file "stats.nc" |
---|
| 1454 | c which can later be used to make the statistic files of the run: |
---|
| 1455 | c "stats") only possible in 3D runs ! |
---|
| 1456 | |
---|
| 1457 | IF (callstats) THEN |
---|
| 1458 | |
---|
[226] | 1459 | call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) |
---|
| 1460 | call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf) |
---|
[520] | 1461 | c call wstats(ngrid,"co2ice","CO2 ice cover", |
---|
| 1462 | c & "kg.m-2",2,co2ice) |
---|
| 1463 | c call wstats(ngrid,"fluxsurf_lw", |
---|
| 1464 | c & "Thermal IR radiative flux to surface","W.m-2",2, |
---|
| 1465 | c & fluxsurf_lw) |
---|
| 1466 | c call wstats(ngrid,"fluxsurf_sw", |
---|
| 1467 | c & "Solar radiative flux to surface","W.m-2",2, |
---|
| 1468 | c & fluxsurf_sw_tot) |
---|
| 1469 | c call wstats(ngrid,"fluxtop_lw", |
---|
| 1470 | c & "Thermal IR radiative flux to space","W.m-2",2, |
---|
| 1471 | c & fluxtop_lw) |
---|
| 1472 | c call wstats(ngrid,"fluxtop_sw", |
---|
| 1473 | c & "Solar radiative flux to space","W.m-2",2, |
---|
| 1474 | c & fluxtop_sw_tot) |
---|
| 1475 | c call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) |
---|
| 1476 | c call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) |
---|
| 1477 | c call wstats(ngrid,"v","Meridional (North-South) wind", |
---|
| 1478 | c & "m.s-1",3,zv) |
---|
| 1479 | c call wstats(ngrid,"w","Vertical (down-up) wind", |
---|
| 1480 | c & "m.s-1",3,pw) |
---|
| 1481 | c call wstats(ngrid,"rho","Atmospheric density","none",3,rho) |
---|
| 1482 | c call wstats(ngrid,"pressure","Pressure","Pa",3,pplay) |
---|
[226] | 1483 | c call wstats(ngrid,"q2", |
---|
| 1484 | c & "Boundary layer eddy kinetic energy", |
---|
| 1485 | c & "m2.s-2",3,q2) |
---|
| 1486 | c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, |
---|
| 1487 | c & emis) |
---|
| 1488 | c call wstats(ngrid,"ssurf","Surface stress","N.m-2", |
---|
| 1489 | c & 2,zstress) |
---|
| 1490 | c call wstats(ngrid,"sw_htrt","sw heat.rate", |
---|
| 1491 | c & "W.m-2",3,zdtsw) |
---|
| 1492 | c call wstats(ngrid,"lw_htrt","lw heat.rate", |
---|
| 1493 | c & "W.m-2",3,zdtlw) |
---|
[42] | 1494 | |
---|
[226] | 1495 | if (tracer) then |
---|
| 1496 | if (water) then |
---|
[520] | 1497 | vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap) |
---|
| 1498 | & *mugaz/mmol(igcm_h2o_vap) |
---|
| 1499 | call wstats(ngrid,"vmr_h2ovapor", |
---|
| 1500 | & "H2O vapor volume mixing ratio","mol/mol", |
---|
| 1501 | & 3,vmr) |
---|
| 1502 | vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) |
---|
| 1503 | & *mugaz/mmol(igcm_h2o_ice) |
---|
| 1504 | call wstats(ngrid,"vmr_h2oice", |
---|
| 1505 | & "H2O ice volume mixing ratio","mol/mol", |
---|
| 1506 | & 3,vmr) |
---|
[226] | 1507 | call wstats(ngrid,"h2o_ice_s", |
---|
| 1508 | & "surface h2o_ice","kg/m2", |
---|
| 1509 | & 2,qsurf(1,igcm_h2o_ice)) |
---|
[520] | 1510 | call wstats(ngrid,'albedo', |
---|
| 1511 | & 'albedo', |
---|
| 1512 | & '',2,albedo(1:ngridmx,1)) |
---|
[226] | 1513 | call wstats(ngrid,"mtot", |
---|
| 1514 | & "total mass of water vapor","kg/m2", |
---|
| 1515 | & 2,mtot) |
---|
| 1516 | call wstats(ngrid,"icetot", |
---|
| 1517 | & "total mass of water ice","kg/m2", |
---|
| 1518 | & 2,icetot) |
---|
[520] | 1519 | call wstats(ngrid,"reffice", |
---|
| 1520 | & "Mean reff","m", |
---|
| 1521 | & 2,rave) |
---|
| 1522 | call wstats(ngrid,"ccntot", |
---|
| 1523 | & "condensation nuclei","Nbr/m2", |
---|
| 1524 | & 2,ccntot) |
---|
| 1525 | call wstats(ngrid,"rice", |
---|
| 1526 | & "Ice particle size","m", |
---|
| 1527 | & 3,rice) |
---|
[226] | 1528 | if (.not.activice) then |
---|
| 1529 | call wstats(ngrid,"tauTESap", |
---|
| 1530 | & "tau abs 825 cm-1","", |
---|
| 1531 | & 2,tauTES) |
---|
[520] | 1532 | else |
---|
| 1533 | call wstats(ngridmx,'tauTES', |
---|
| 1534 | & 'tau abs 825 cm-1', |
---|
| 1535 | & '',2,taucloudtes) |
---|
[226] | 1536 | endif |
---|
| 1537 | |
---|
| 1538 | endif ! of if (water) |
---|
| 1539 | |
---|
| 1540 | if (thermochem.or.photochem) then |
---|
| 1541 | do iq=1,nq |
---|
[520] | 1542 | if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or. |
---|
| 1543 | . (noms(iq).eq."co").or.(noms(iq).eq."n2").or. |
---|
| 1544 | . (noms(iq).eq."h2").or. |
---|
| 1545 | . (noms(iq).eq."o3")) then |
---|
| 1546 | do l=1,nlayer |
---|
| 1547 | do ig=1,ngrid |
---|
| 1548 | vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq) |
---|
| 1549 | end do |
---|
| 1550 | end do |
---|
| 1551 | call wstats(ngrid,"vmr_"//trim(noms(iq)), |
---|
| 1552 | . "Volume mixing ratio","mol/mol",3,vmr) |
---|
| 1553 | endif |
---|
| 1554 | ! do ig = 1,ngrid |
---|
| 1555 | ! colden(ig,iq) = 0. !FL |
---|
| 1556 | ! end do |
---|
| 1557 | ! do l=1,nlayer !FL |
---|
| 1558 | ! do ig=1,ngrid !FL |
---|
| 1559 | ! colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq) !FL |
---|
| 1560 | ! $ *(pplev(ig,l)-pplev(ig,l+1)) !FL |
---|
| 1561 | ! $ *6.022e22/(mmol(iq)*g) !FL |
---|
| 1562 | ! end do !FL |
---|
| 1563 | ! end do !FL |
---|
| 1564 | ! call wstats(ngrid,"c_"//trim(noms(iq)), !FL |
---|
| 1565 | ! $ "column","mol cm-2",2,colden(1,iq)) !FL |
---|
[334] | 1566 | end do |
---|
| 1567 | end if ! of if (thermochem.or.photochem) |
---|
[226] | 1568 | |
---|
[334] | 1569 | end if ! of if (tracer) |
---|
[226] | 1570 | |
---|
| 1571 | IF(lastcall) THEN |
---|
| 1572 | write (*,*) "Writing stats..." |
---|
| 1573 | call mkstats(ierr) |
---|
| 1574 | ENDIF |
---|
| 1575 | |
---|
[42] | 1576 | ENDIF !if callstats |
---|
| 1577 | |
---|
| 1578 | c (Store EOF for Mars Climate database software) |
---|
| 1579 | IF (calleofdump) THEN |
---|
| 1580 | CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps) |
---|
| 1581 | ENDIF |
---|
| 1582 | |
---|
[234] | 1583 | |
---|
[226] | 1584 | #ifdef MESOSCALE |
---|
[234] | 1585 | !!! |
---|
| 1586 | !!! OUTPUT FIELDS |
---|
| 1587 | !!! |
---|
[81] | 1588 | wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature |
---|
| 1589 | wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice |
---|
[82] | 1590 | mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice |
---|
[308] | 1591 | !! JF |
---|
| 1592 | TAU_lay(:)=tau(:,1)!!true opacity (not a reference like tauref) |
---|
| 1593 | IF (igcm_dust_mass .ne. 0) THEN |
---|
| 1594 | qsurfice_dust(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass) |
---|
| 1595 | ENDIF |
---|
[81] | 1596 | IF (igcm_h2o_ice .ne. 0) THEN |
---|
| 1597 | qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice) |
---|
| 1598 | vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
| 1599 | . * mugaz / mmol(igcm_h2o_ice) |
---|
| 1600 | ENDIF |
---|
[308] | 1601 | !! Dust quantity integration along the vertical axe |
---|
| 1602 | dustot(:)=0 |
---|
| 1603 | do ig=1,ngrid |
---|
| 1604 | do l=1,nlayermx |
---|
| 1605 | dustot(ig) = dustot(ig) + |
---|
| 1606 | & zq(ig,l,igcm_dust_mass) |
---|
| 1607 | & * (pplev(ig,l) - pplev(ig,l+1)) / g |
---|
| 1608 | enddo |
---|
| 1609 | enddo |
---|
[226] | 1610 | c AUTOMATICALLY GENERATED FROM REGISTRY |
---|
[42] | 1611 | #include "fill_save.inc" |
---|
[226] | 1612 | #else |
---|
[528] | 1613 | #ifndef MESOINI |
---|
[42] | 1614 | |
---|
[226] | 1615 | c ========================================================== |
---|
| 1616 | c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing |
---|
| 1617 | c any variable for diagnostic (output with period |
---|
| 1618 | c "ecritphy", set in "run.def") |
---|
| 1619 | c ========================================================== |
---|
| 1620 | c WRITEDIAGFI can ALSO be called from any other subroutines |
---|
| 1621 | c for any variables !! |
---|
| 1622 | c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, |
---|
| 1623 | c & emis) |
---|
[520] | 1624 | c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay) |
---|
| 1625 | c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) |
---|
[226] | 1626 | call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, |
---|
| 1627 | & tsurf) |
---|
| 1628 | call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) |
---|
| 1629 | call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, |
---|
| 1630 | & co2ice) |
---|
[284] | 1631 | |
---|
[226] | 1632 | c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7", |
---|
| 1633 | c & "K",2,zt(1,7)) |
---|
[284] | 1634 | c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, |
---|
| 1635 | c & fluxsurf_lw) |
---|
| 1636 | c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2, |
---|
| 1637 | c & fluxsurf_sw_tot) |
---|
| 1638 | c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2, |
---|
| 1639 | c & fluxtop_lw) |
---|
| 1640 | c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2, |
---|
| 1641 | c & fluxtop_sw_tot) |
---|
[520] | 1642 | c call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) |
---|
| 1643 | c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) |
---|
| 1644 | c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) |
---|
| 1645 | c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) |
---|
| 1646 | c call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) |
---|
[226] | 1647 | c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2) |
---|
[520] | 1648 | c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) |
---|
[226] | 1649 | c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay) |
---|
| 1650 | c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, |
---|
| 1651 | c & zstress) |
---|
| 1652 | c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate', |
---|
| 1653 | c & 'w.m-2',3,zdtsw) |
---|
| 1654 | c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', |
---|
| 1655 | c & 'w.m-2',3,zdtlw) |
---|
[520] | 1656 | if (.not.activice) then |
---|
| 1657 | CALL WRITEDIAGFI(ngridmx,'tauTESap', |
---|
| 1658 | & 'tau abs 825 cm-1', |
---|
| 1659 | & '',2,tauTES) |
---|
| 1660 | else |
---|
| 1661 | CALL WRITEDIAGFI(ngridmx,'tauTES', |
---|
| 1662 | & 'tau abs 825 cm-1', |
---|
| 1663 | & '',2,taucloudtes) |
---|
| 1664 | endif |
---|
[42] | 1665 | |
---|
[528] | 1666 | #else |
---|
| 1667 | !!! this is to ensure correct initialisation of mesoscale model |
---|
| 1668 | call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, |
---|
| 1669 | & tsurf) |
---|
| 1670 | call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) |
---|
| 1671 | call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, |
---|
| 1672 | & co2ice) |
---|
| 1673 | call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) |
---|
| 1674 | call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) |
---|
| 1675 | call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) |
---|
[299] | 1676 | call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, |
---|
| 1677 | & emis) |
---|
| 1678 | call WRITEDIAGFI(ngrid,"tsoil","Soil temperature", |
---|
| 1679 | & "K",3,tsoil) |
---|
| 1680 | call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia", |
---|
| 1681 | & "K",3,inertiedat) |
---|
| 1682 | #endif |
---|
[42] | 1683 | |
---|
[299] | 1684 | |
---|
[226] | 1685 | c ---------------------------------------------------------- |
---|
| 1686 | c Outputs of the CO2 cycle |
---|
| 1687 | c ---------------------------------------------------------- |
---|
[42] | 1688 | |
---|
[226] | 1689 | if (tracer.and.(igcm_co2.ne.0)) then |
---|
| 1690 | ! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer", |
---|
| 1691 | ! & "kg/kg",2,zq(1,1,igcm_co2)) |
---|
| 1692 | ! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio", |
---|
| 1693 | ! & "kg/kg",3,zq(1,1,igcm_co2)) |
---|
| 1694 | |
---|
| 1695 | ! Compute co2 column |
---|
[286] | 1696 | co2col(:)=0 |
---|
[226] | 1697 | do l=1,nlayermx |
---|
| 1698 | do ig=1,ngrid |
---|
| 1699 | co2col(ig)=co2col(ig)+ |
---|
| 1700 | & zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g |
---|
| 1701 | enddo |
---|
| 1702 | enddo |
---|
| 1703 | call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2, |
---|
| 1704 | & co2col) |
---|
[520] | 1705 | !!!!! FL |
---|
| 1706 | ! do iq = 1,nq |
---|
| 1707 | ! if (noms(iq) .ne. "dust_mass" .and. |
---|
| 1708 | ! $ noms(iq) .ne. "dust_number") then |
---|
| 1709 | ! call writediagfi(ngrid,"c_"//trim(noms(iq)), |
---|
| 1710 | ! $ "column","mol cm-2",2,colden(1,iq)) |
---|
| 1711 | ! end if |
---|
| 1712 | ! end do |
---|
| 1713 | !!!!! FL |
---|
[226] | 1714 | endif ! of if (tracer.and.(igcm_co2.ne.0)) |
---|
[42] | 1715 | |
---|
[226] | 1716 | c ---------------------------------------------------------- |
---|
| 1717 | c Outputs of the water cycle |
---|
| 1718 | c ---------------------------------------------------------- |
---|
| 1719 | if (tracer) then |
---|
| 1720 | if (water) then |
---|
| 1721 | |
---|
[299] | 1722 | #ifdef MESOINI |
---|
| 1723 | !!!! waterice = q01, voir readmeteo.F90 |
---|
| 1724 | call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice), |
---|
| 1725 | & 'kg/kg',3, |
---|
| 1726 | & zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)) |
---|
| 1727 | !!!! watervapor = q02, voir readmeteo.F90 |
---|
| 1728 | call WRITEDIAGFI(ngridmx,'q02',noms(igcm_h2o_vap), |
---|
| 1729 | & 'kg/kg',3, |
---|
| 1730 | & zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)) |
---|
| 1731 | !!!! surface waterice qsurf02 (voir readmeteo) |
---|
| 1732 | call WRITEDIAGFI(ngridmx,'qsurf02','surface tracer', |
---|
| 1733 | & 'kg.m-2',2, |
---|
| 1734 | & qsurf(1:ngridmx,igcm_h2o_ice)) |
---|
| 1735 | #endif |
---|
| 1736 | |
---|
[226] | 1737 | CALL WRITEDIAGFI(ngridmx,'mtot', |
---|
| 1738 | & 'total mass of water vapor', |
---|
| 1739 | & 'kg/m2',2,mtot) |
---|
| 1740 | CALL WRITEDIAGFI(ngridmx,'icetot', |
---|
| 1741 | & 'total mass of water ice', |
---|
| 1742 | & 'kg/m2',2,icetot) |
---|
[411] | 1743 | c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) |
---|
| 1744 | c & *mugaz/mmol(igcm_h2o_ice) |
---|
| 1745 | c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr', |
---|
| 1746 | c & 'mol/mol',3,vmr) |
---|
| 1747 | c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap) |
---|
| 1748 | c & *mugaz/mmol(igcm_h2o_vap) |
---|
| 1749 | c call WRITEDIAGFI(ngridmx,'vmr_h2ovap','h2o vap vmr', |
---|
| 1750 | c & 'mol/mol',3,vmr) |
---|
| 1751 | c CALL WRITEDIAGFI(ngridmx,'reffice', |
---|
| 1752 | c & 'Mean reff', |
---|
| 1753 | c & 'm',2,rave) |
---|
[520] | 1754 | c CALL WRITEDIAGFI(ngrid,"ccntot", |
---|
| 1755 | c & "condensation nuclei","Nbr/m2", |
---|
| 1756 | c & 2,ccntot) |
---|
[411] | 1757 | c call WRITEDIAGFI(ngridmx,'rice','Ice particle size', |
---|
| 1758 | c & 'm',3,rice) |
---|
[226] | 1759 | call WRITEDIAGFI(ngridmx,'h2o_ice_s', |
---|
| 1760 | & 'surface h2o_ice', |
---|
| 1761 | & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) |
---|
[283] | 1762 | |
---|
[520] | 1763 | c if (caps) then |
---|
| 1764 | c do ig=1,ngridmx |
---|
| 1765 | c if (watercaptag(ig)) watercapflag(ig) = 1 |
---|
| 1766 | c enddo |
---|
[358] | 1767 | c CALL WRITEDIAGFI(ngridmx,'watercaptag', |
---|
[520] | 1768 | c & 'Ice water caps', |
---|
| 1769 | c & '',2,watercapflag) |
---|
| 1770 | cs endif |
---|
[358] | 1771 | c CALL WRITEDIAGFI(ngridmx,'albedo', |
---|
| 1772 | c & 'albedo', |
---|
| 1773 | c & '',2,albedo(1:ngridmx,1)) |
---|
[226] | 1774 | endif !(water) |
---|
| 1775 | |
---|
| 1776 | |
---|
| 1777 | if (water.and..not.photochem) then |
---|
| 1778 | iq=nq |
---|
| 1779 | c write(str2(1:2),'(i2.2)') iq |
---|
| 1780 | c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud', |
---|
| 1781 | c & 'kg.m-2',2,zdqscloud(1,iq)) |
---|
| 1782 | c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim', |
---|
| 1783 | c & 'kg/kg',3,zdqchim(1,1,iq)) |
---|
| 1784 | c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif', |
---|
| 1785 | c & 'kg/kg',3,zdqdif(1,1,iq)) |
---|
| 1786 | c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj', |
---|
| 1787 | c & 'kg/kg',3,zdqadj(1,1,iq)) |
---|
| 1788 | c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c', |
---|
| 1789 | c & 'kg/kg',3,zdqc(1,1,iq)) |
---|
| 1790 | endif !(water.and..not.photochem) |
---|
| 1791 | endif |
---|
| 1792 | |
---|
| 1793 | c ---------------------------------------------------------- |
---|
| 1794 | c Outputs of the dust cycle |
---|
| 1795 | c ---------------------------------------------------------- |
---|
| 1796 | |
---|
| 1797 | c call WRITEDIAGFI(ngridmx,'tauref', |
---|
| 1798 | c & 'Dust ref opt depth','NU',2,tauref) |
---|
| 1799 | |
---|
| 1800 | if (tracer.and.(dustbin.ne.0)) then |
---|
| 1801 | c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1)) |
---|
| 1802 | if (doubleq) then |
---|
[411] | 1803 | c call WRITEDIAGFI(ngridmx,'qsurf','qsurf', |
---|
| 1804 | c & 'kg.m-2',2,qsurf(1,igcm_dust_mass)) |
---|
| 1805 | c call WRITEDIAGFI(ngridmx,'Nsurf','N particles', |
---|
| 1806 | c & 'N.m-2',2,qsurf(1,igcm_dust_number)) |
---|
[226] | 1807 | c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift', |
---|
| 1808 | c & 'kg.m-2.s-1',2,zdqsdev(1,1)) |
---|
[411] | 1809 | c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation', |
---|
| 1810 | c & 'kg.m-2.s-1',2,zdqssed(1,1)) |
---|
| 1811 | c call WRITEDIAGFI(ngridmx,'dqsdif','diffusion', |
---|
| 1812 | c & 'kg.m-2.s-1',2,zdqsdif(1,1)) |
---|
| 1813 | c call WRITEDIAGFI(ngridmx,'reffdust','reffdust', |
---|
| 1814 | c & 'm',3,rdust*ref_r0) |
---|
| 1815 | c call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr', |
---|
| 1816 | c & 'kg/kg',3,pq(1,1,igcm_dust_mass)) |
---|
| 1817 | c call WRITEDIAGFI(ngridmx,'dustN','Dust number', |
---|
| 1818 | c & 'part/kg',3,pq(1,1,igcm_dust_number)) |
---|
[299] | 1819 | #ifdef MESOINI |
---|
| 1820 | call WRITEDIAGFI(ngridmx,'dustN','Dust number', |
---|
[520] | 1821 | & 'part/kg',3,pq(1,1,igcm_dust_number)) |
---|
[299] | 1822 | #endif |
---|
[226] | 1823 | else |
---|
| 1824 | do iq=1,dustbin |
---|
| 1825 | write(str2(1:2),'(i2.2)') iq |
---|
| 1826 | call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio', |
---|
| 1827 | & 'kg/kg',3,zq(1,1,iq)) |
---|
| 1828 | call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf', |
---|
| 1829 | & 'kg.m-2',2,qsurf(1,iq)) |
---|
| 1830 | end do |
---|
| 1831 | endif ! (doubleq) |
---|
[358] | 1832 | |
---|
| 1833 | if (scavenging) then |
---|
[520] | 1834 | c call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr', |
---|
| 1835 | c & 'kg/kg',3,pq(1,1,igcm_ccn_mass)) |
---|
| 1836 | c call WRITEDIAGFI(ngridmx,'ccnN','CCN number', |
---|
| 1837 | c & 'part/kg',3,pq(1,1,igcm_ccn_number)) |
---|
[358] | 1838 | endif ! (scavenging) |
---|
| 1839 | |
---|
[226] | 1840 | c if (submicron) then |
---|
| 1841 | c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr', |
---|
| 1842 | c & 'kg/kg',3,pq(1,1,igcm_dust_submicron)) |
---|
| 1843 | c endif ! (submicron) |
---|
| 1844 | end if ! (tracer.and.(dustbin.ne.0)) |
---|
| 1845 | |
---|
| 1846 | c ---------------------------------------------------------- |
---|
[267] | 1847 | c ---------------------------------------------------------- |
---|
| 1848 | c PBL OUTPUS |
---|
| 1849 | c ---------------------------------------------------------- |
---|
| 1850 | c ---------------------------------------------------------- |
---|
| 1851 | |
---|
| 1852 | |
---|
| 1853 | c ---------------------------------------------------------- |
---|
| 1854 | c Outputs of surface layer |
---|
| 1855 | c ---------------------------------------------------------- |
---|
| 1856 | |
---|
| 1857 | |
---|
[353] | 1858 | z_out=0. |
---|
[267] | 1859 | if (calltherm .and. (z_out .gt. 0.)) then |
---|
| 1860 | |
---|
[529] | 1861 | call pbl_parameters(ngrid,nlayer,ps,zplay,z0, |
---|
| 1862 | & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, |
---|
[528] | 1863 | & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) |
---|
[267] | 1864 | |
---|
| 1865 | else |
---|
| 1866 | if((.not. calltherm).and.(z_out .gt. 0.)) then |
---|
| 1867 | print*, 'WARNING : no interpolation in surface-layer :' |
---|
| 1868 | print*, 'Outputing surface-layer quantities without thermals |
---|
| 1869 | & does not make sense' |
---|
| 1870 | endif |
---|
| 1871 | endif |
---|
| 1872 | |
---|
| 1873 | c ---------------------------------------------------------- |
---|
[226] | 1874 | c Outputs of thermals |
---|
| 1875 | c ---------------------------------------------------------- |
---|
| 1876 | if (calltherm) then |
---|
| 1877 | |
---|
| 1878 | ! call WRITEDIAGFI(ngrid,'dtke', |
---|
| 1879 | ! & 'tendance tke thermiques','m**2/s**2', |
---|
| 1880 | ! & 3,dtke_th) |
---|
| 1881 | ! call WRITEDIAGFI(ngrid,'d_u_ajs', |
---|
| 1882 | ! & 'tendance u thermiques','m/s', |
---|
| 1883 | ! & 3,pdu_th*ptimestep) |
---|
| 1884 | ! call WRITEDIAGFI(ngrid,'d_v_ajs', |
---|
| 1885 | ! & 'tendance v thermiques','m/s', |
---|
| 1886 | ! & 3,pdv_th*ptimestep) |
---|
| 1887 | ! if (tracer) then |
---|
| 1888 | ! if (nq .eq. 2) then |
---|
| 1889 | ! call WRITEDIAGFI(ngrid,'deltaq_th', |
---|
| 1890 | ! & 'delta q thermiques','kg/kg', |
---|
| 1891 | ! & 3,ptimestep*pdq_th(:,:,2)) |
---|
| 1892 | ! endif |
---|
| 1893 | ! endif |
---|
| 1894 | |
---|
[523] | 1895 | call WRITEDIAGFI(ngridmx,'zmax_th', |
---|
| 1896 | & 'hauteur du thermique','m', |
---|
| 1897 | & 2,zmax_th) |
---|
[226] | 1898 | call WRITEDIAGFI(ngridmx,'hfmax_th', |
---|
| 1899 | & 'maximum TH heat flux','K.m/s', |
---|
| 1900 | & 2,hfmax_th) |
---|
[528] | 1901 | call WRITEDIAGFI(ngridmx,'wstar', |
---|
[226] | 1902 | & 'maximum TH vertical velocity','m/s', |
---|
[528] | 1903 | & 2,wstar) |
---|
[226] | 1904 | |
---|
| 1905 | endif |
---|
| 1906 | |
---|
| 1907 | c ---------------------------------------------------------- |
---|
[267] | 1908 | c ---------------------------------------------------------- |
---|
| 1909 | c END OF PBL OUTPUS |
---|
| 1910 | c ---------------------------------------------------------- |
---|
| 1911 | c ---------------------------------------------------------- |
---|
| 1912 | |
---|
| 1913 | |
---|
| 1914 | c ---------------------------------------------------------- |
---|
[226] | 1915 | c Output in netcdf file "diagsoil.nc" for subterranean |
---|
| 1916 | c variables (output every "ecritphy", as for writediagfi) |
---|
| 1917 | c ---------------------------------------------------------- |
---|
| 1918 | |
---|
| 1919 | ! Write soil temperature |
---|
| 1920 | ! call writediagsoil(ngrid,"soiltemp","Soil temperature","K", |
---|
| 1921 | ! & 3,tsoil) |
---|
| 1922 | ! Write surface temperature |
---|
| 1923 | ! call writediagsoil(ngrid,"tsurf","Surface temperature","K", |
---|
| 1924 | ! & 2,tsurf) |
---|
| 1925 | |
---|
| 1926 | c ========================================================== |
---|
| 1927 | c END OF WRITEDIAGFI |
---|
| 1928 | c ========================================================== |
---|
| 1929 | #endif |
---|
| 1930 | |
---|
[42] | 1931 | ELSE ! if(ngrid.eq.1) |
---|
| 1932 | |
---|
| 1933 | print*,'Ls =',zls*180./pi, |
---|
| 1934 | & ' tauref(700 Pa) =',tauref |
---|
| 1935 | c ---------------------------------------------------------------------- |
---|
| 1936 | c Output in grads file "g1d" (ONLY when using testphys1d) |
---|
| 1937 | c (output at every X physical timestep) |
---|
| 1938 | c ---------------------------------------------------------------------- |
---|
| 1939 | c |
---|
| 1940 | c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2') |
---|
| 1941 | c CALL writeg1d(ngrid,1,tsurf,'tsurf','K') |
---|
| 1942 | c CALL writeg1d(ngrid,1,ps,'ps','Pa') |
---|
| 1943 | |
---|
| 1944 | c CALL writeg1d(ngrid,nlayer,zt,'T','K') |
---|
| 1945 | c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1') |
---|
| 1946 | c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1') |
---|
| 1947 | c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1') |
---|
| 1948 | |
---|
[226] | 1949 | ! THERMALS STUFF 1D |
---|
[42] | 1950 | |
---|
[529] | 1951 | z_out=1. |
---|
[268] | 1952 | if (calltherm .and. (z_out .gt. 0.)) then |
---|
| 1953 | |
---|
[529] | 1954 | call pbl_parameters(ngrid,nlayer,ps,zplay,z0, |
---|
| 1955 | & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, |
---|
[528] | 1956 | & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) |
---|
[268] | 1957 | |
---|
| 1958 | else |
---|
| 1959 | if((.not. calltherm).and.(z_out .gt. 0.)) then |
---|
| 1960 | print*, 'WARNING : no interpolation in surface-layer :' |
---|
| 1961 | print*, 'Outputing surface-layer quantities without thermals |
---|
| 1962 | & does not make sense' |
---|
| 1963 | endif |
---|
| 1964 | endif |
---|
| 1965 | |
---|
[226] | 1966 | if(calltherm) then |
---|
| 1967 | |
---|
| 1968 | call WRITEDIAGFI(ngridmx,'lmax_th', |
---|
| 1969 | & 'hauteur du thermique','point', |
---|
| 1970 | & 0,lmax_th_out) |
---|
[528] | 1971 | call WRITEDIAGFI(ngridmx,'zmax_th', |
---|
| 1972 | & 'hauteur du thermique','m', |
---|
| 1973 | & 0,zmax_th) |
---|
[226] | 1974 | call WRITEDIAGFI(ngridmx,'hfmax_th', |
---|
| 1975 | & 'maximum TH heat flux','K.m/s', |
---|
| 1976 | & 0,hfmax_th) |
---|
[528] | 1977 | call WRITEDIAGFI(ngridmx,'wstar', |
---|
[226] | 1978 | & 'maximum TH vertical velocity','m/s', |
---|
[528] | 1979 | & 0,wstar) |
---|
[226] | 1980 | |
---|
| 1981 | co2col(:)=0. |
---|
| 1982 | if (tracer) then |
---|
| 1983 | do l=1,nlayermx |
---|
| 1984 | do ig=1,ngrid |
---|
| 1985 | co2col(ig)=co2col(ig)+ |
---|
| 1986 | & zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g |
---|
| 1987 | enddo |
---|
| 1988 | enddo |
---|
| 1989 | |
---|
| 1990 | end if |
---|
| 1991 | call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' & |
---|
| 1992 | & ,'kg/m-2',0,co2col) |
---|
| 1993 | endif |
---|
| 1994 | call WRITEDIAGFI(ngrid,'w','vertical velocity' & |
---|
| 1995 | & ,'m/s',1,pw) |
---|
| 1996 | call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2) |
---|
| 1997 | call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0, |
---|
| 1998 | & tsurf) |
---|
[277] | 1999 | call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu) |
---|
| 2000 | call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv) |
---|
[226] | 2001 | |
---|
| 2002 | call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay) |
---|
| 2003 | call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev) |
---|
[358] | 2004 | call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho) |
---|
[544] | 2005 | ! call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", & |
---|
| 2006 | ! & "K.s-1",1,dtrad/zpopsk) |
---|
| 2007 | ! call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate', |
---|
| 2008 | ! & 'w.m-2',1,zdtsw/zpopsk) |
---|
| 2009 | ! call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', |
---|
| 2010 | ! & 'w.m-2',1,zdtlw/zpopsk) |
---|
| 2011 | |
---|
[226] | 2012 | ! or output in diagfi.nc (for testphys1d) |
---|
| 2013 | call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps) |
---|
| 2014 | call WRITEDIAGFI(ngridmx,'temp','Temperature', |
---|
| 2015 | & 'K',1,zt) |
---|
[358] | 2016 | |
---|
[226] | 2017 | if(tracer) then |
---|
| 2018 | c CALL writeg1d(ngrid,1,tau,'tau','SI') |
---|
| 2019 | do iq=1,nq |
---|
| 2020 | c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') |
---|
| 2021 | call WRITEDIAGFI(ngridmx,trim(noms(iq)), |
---|
| 2022 | & trim(noms(iq)),'kg/kg',1,zq(1,1,iq)) |
---|
| 2023 | end do |
---|
[358] | 2024 | if (doubleq) then |
---|
| 2025 | call WRITEDIAGFI(ngridmx,'rdust','rdust', |
---|
| 2026 | & 'm',1,rdust) |
---|
| 2027 | endif |
---|
[226] | 2028 | end if |
---|
[358] | 2029 | |
---|
[520] | 2030 | cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc |
---|
[358] | 2031 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
[520] | 2032 | IF (water) THEN |
---|
| 2033 | CALL WRITEDIAGFI(ngridmx,'tauTESap', |
---|
[358] | 2034 | & 'tau abs 825 cm-1', |
---|
[520] | 2035 | & '',0,tauTES) |
---|
[358] | 2036 | |
---|
[520] | 2037 | CALL WRITEDIAGFI(ngridmx,'tauTES', |
---|
| 2038 | & 'tau abs 825 cm-1', |
---|
| 2039 | & '',0,taucloudtes) |
---|
| 2040 | |
---|
| 2041 | mtot = 0 |
---|
| 2042 | icetot = 0 |
---|
| 2043 | h2otot = qsurf(1,igcm_h2o_ice) |
---|
| 2044 | rave = 0 |
---|
[358] | 2045 | do l=1,nlayer |
---|
[520] | 2046 | mtot = mtot + zq(1,l,igcm_h2o_vap) |
---|
| 2047 | & * (pplev(1,l) - pplev(1,l+1)) / g |
---|
| 2048 | icetot = icetot + zq(1,l,igcm_h2o_ice) |
---|
| 2049 | & * (pplev(1,l) - pplev(1,l+1)) / g |
---|
| 2050 | rave = rave + zq(1,l,igcm_h2o_ice) |
---|
| 2051 | & * (pplev(1,l) - pplev(1,l+1)) / g |
---|
| 2052 | & * rice(1,l) * (1.+nuice_ref) |
---|
[358] | 2053 | end do |
---|
[520] | 2054 | rave=rave/max(icetot,1.e-30) |
---|
| 2055 | h2otot = h2otot+mtot+icetot |
---|
| 2056 | |
---|
| 2057 | |
---|
| 2058 | if (scavenging) then |
---|
| 2059 | ccntot= 0 |
---|
| 2060 | call watersat(ngridmx*nlayermx,zt,pplay,zqsat) |
---|
| 2061 | do l=1,nlayermx |
---|
| 2062 | ccntot = ccntot + |
---|
| 2063 | & zq(1,l,igcm_ccn_number)*tauscaling(1) |
---|
| 2064 | & *(pplev(1,l) - pplev(1,l+1)) / g |
---|
| 2065 | satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l) |
---|
| 2066 | satu(1,l) = (max(satu(1,l),float(1))-1) |
---|
| 2067 | ! & * zq(1,l,igcm_h2o_vap) * |
---|
| 2068 | ! & (pplev(1,l) - pplev(1,l+1)) / g |
---|
| 2069 | enddo |
---|
| 2070 | |
---|
| 2071 | CALL WRITEDIAGFI(ngridmx,'ccntot', |
---|
| 2072 | & 'ccntot', |
---|
| 2073 | & 'nbr/m2',0,ccntot) |
---|
| 2074 | endif |
---|
| 2075 | |
---|
| 2076 | |
---|
| 2077 | CALL WRITEDIAGFI(ngridmx,'h2otot', |
---|
| 2078 | & 'h2otot', |
---|
| 2079 | & 'kg/m2',0,h2otot) |
---|
| 2080 | CALL WRITEDIAGFI(ngridmx,'mtot', |
---|
| 2081 | & 'mtot', |
---|
| 2082 | & 'kg/m2',0,mtot) |
---|
| 2083 | CALL WRITEDIAGFI(ngridmx,'icetot', |
---|
| 2084 | & 'icetot', |
---|
| 2085 | & 'kg/m2',0,icetot) |
---|
| 2086 | CALL WRITEDIAGFI(ngridmx,'reffice', |
---|
| 2087 | & 'reffice', |
---|
| 2088 | & 'm',0,rave) |
---|
[411] | 2089 | |
---|
| 2090 | |
---|
| 2091 | do iq=1,nq |
---|
| 2092 | call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_s', |
---|
| 2093 | & trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq)) |
---|
| 2094 | end do |
---|
| 2095 | |
---|
[358] | 2096 | |
---|
[520] | 2097 | call WRITEDIAGFI(ngridmx,'zdqsed_dustq','sedimentation q', |
---|
| 2098 | & 'kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass)) |
---|
| 2099 | call WRITEDIAGFI(ngridmx,'zdqsed_dustN','sedimentation N', |
---|
| 2100 | & 'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_dust_number)) |
---|
[358] | 2101 | |
---|
[520] | 2102 | call WRITEDIAGFI(ngridmx,'zdqcloud*dt ice','cloud ice', |
---|
| 2103 | & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)*ptimestep) |
---|
| 2104 | call WRITEDIAGFI(ngridmx,'zdqcloud*dt vap','cloud vap', |
---|
| 2105 | & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap)*ptimestep) |
---|
| 2106 | call WRITEDIAGFI(ngridmx,'zdqdif*dt ice','dif ice', |
---|
| 2107 | & 'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_ice)*ptimestep) |
---|
| 2108 | call WRITEDIAGFI(ngridmx,'zdqdif*dt vap','dif vap', |
---|
| 2109 | & 'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_vap)*ptimestep) |
---|
| 2110 | call WRITEDIAGFI(ngridmx,'zdqdif*dt vap 0','dif vap', |
---|
| 2111 | & 'kg.m-2.s-1',0,zdqdif(1,1,igcm_h2o_vap)*ptimestep) |
---|
| 2112 | |
---|
| 2113 | ! if (scavenging) then |
---|
| 2114 | ! call WRITEDIAGFI(ngridmx,'zdqsed_ccnq','sedimentation q', |
---|
| 2115 | ! & 'kg.m-2.s-1',1,zdqsed(1,:,igcm_ccn_mass)) |
---|
| 2116 | ! call WRITEDIAGFI(ngridmx,'zdqsed_ccnN','sedimentation N', |
---|
| 2117 | ! & 'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_ccn_number)) |
---|
| 2118 | ! endif |
---|
| 2119 | ! call WRITEDIAGFI(ngridmx,'zdqsed_ice','sedimentation q', |
---|
| 2120 | ! & 'kg.m-2.s-1',1,zdqsed(1,:,igcm_h2o_ice)) |
---|
| 2121 | ! |
---|
[411] | 2122 | |
---|
[520] | 2123 | call WRITEDIAGFI(ngrid,"rice","ice radius","m",1, |
---|
[411] | 2124 | & rice) |
---|
[520] | 2125 | call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1, |
---|
| 2126 | & satu) |
---|
| 2127 | ENDIF ! of IF (water) |
---|
[411] | 2128 | |
---|
[358] | 2129 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
| 2130 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
| 2131 | |
---|
| 2132 | |
---|
[226] | 2133 | zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g |
---|
| 2134 | |
---|
| 2135 | do l=2,nlayer-1 |
---|
| 2136 | tmean=zt(1,l) |
---|
| 2137 | if(zt(1,l).ne.zt(1,l-1)) |
---|
| 2138 | & tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1)) |
---|
| 2139 | zlocal(l)= zlocal(l-1) |
---|
| 2140 | & -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g |
---|
| 2141 | enddo |
---|
| 2142 | zlocal(nlayer)= zlocal(nlayer-1)- |
---|
| 2143 | & log(pplay(1,nlayer)/pplay(1,nlayer-1))* |
---|
| 2144 | & rnew(1,nlayer)*tmean/g |
---|
| 2145 | |
---|
[42] | 2146 | END IF ! if(ngrid.ne.1) |
---|
| 2147 | |
---|
| 2148 | icount=icount+1 |
---|
| 2149 | RETURN |
---|
| 2150 | END |
---|