[1549] | 1 | module physiq_mod |
---|
| 2 | |
---|
| 3 | implicit none |
---|
| 4 | |
---|
| 5 | contains |
---|
| 6 | |
---|
[751] | 7 | subroutine physiq(ngrid,nlayer,nq, & |
---|
[787] | 8 | nametrac, & |
---|
[253] | 9 | firstcall,lastcall, & |
---|
| 10 | pday,ptime,ptimestep, & |
---|
| 11 | pplev,pplay,pphi, & |
---|
| 12 | pu,pv,pt,pq, & |
---|
[1312] | 13 | flxw, & |
---|
[1576] | 14 | pdu,pdv,pdt,pdq,pdpsrf) |
---|
[253] | 15 | |
---|
[726] | 16 | use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind |
---|
[1016] | 17 | use watercommon_h, only : RLVTT, Psat_water,epsi |
---|
[1216] | 18 | use gases_h, only: gnom, gfrac |
---|
[1482] | 19 | use radcommon_h, only: sigma, glat, grav, BWNV |
---|
[1308] | 20 | use radii_mod, only: h2o_reffrad, co2_reffrad |
---|
[1216] | 21 | use aerosol_mod, only: iaero_co2, iaero_h2o |
---|
[1482] | 22 | use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, & |
---|
[1216] | 23 | dryness, watercaptag |
---|
| 24 | use comdiurn_h, only: coslat, sinlat, coslon, sinlon |
---|
[1327] | 25 | use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen |
---|
[1216] | 26 | use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat |
---|
[1543] | 27 | use geometry_mod, only: latitude, longitude, cell_area |
---|
[1542] | 28 | USE comgeomfi_h, only: totarea, totarea_planet |
---|
[1216] | 29 | USE tracer_h, only: noms, mmol, radius, rho_q, qext, & |
---|
| 30 | alpha_lift, alpha_devil, qextrhor, & |
---|
| 31 | igcm_h2o_ice, igcm_h2o_vap, igcm_dustbin, & |
---|
| 32 | igcm_co2_ice |
---|
[1525] | 33 | use time_phylmdz_mod, only: ecritphy, iphysiq, nday |
---|
[1669] | 34 | use phyetat0_mod, only: phyetat0 |
---|
[1216] | 35 | use phyredem, only: physdem0, physdem1 |
---|
[1529] | 36 | use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, & |
---|
| 37 | noceanmx |
---|
[1297] | 38 | use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, & |
---|
[1529] | 39 | ini_surf_heat_transp_mod, & |
---|
[1297] | 40 | ocean_slab_get_vars,ocean_slab_final |
---|
| 41 | use surf_heat_transp_mod,only :init_masquv |
---|
[1295] | 42 | use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval |
---|
| 43 | use mod_phys_lmdz_para, only : is_master |
---|
[1308] | 44 | use planete_mod, only: apoastr, periastr, year_day, peri_day, & |
---|
| 45 | obliquit, nres, z0 |
---|
[1524] | 46 | use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp |
---|
| 47 | use time_phylmdz_mod, only: daysec |
---|
[1397] | 48 | use callkeys_mod |
---|
[1801] | 49 | use conc_mod |
---|
[1836] | 50 | use phys_state_var_mod |
---|
| 51 | use turb_mod, only : q2,sensibFlux,turb_resolved |
---|
| 52 | #ifndef MESOSCALE |
---|
[1622] | 53 | use vertical_layers_mod, only: presnivs, pseudoalt |
---|
[1682] | 54 | use mod_phys_lmdz_omp_data, ONLY: is_omp_master |
---|
[1836] | 55 | #else |
---|
| 56 | use comm_wrf, only : comm_HR_SW, comm_HR_LW, & |
---|
| 57 | comm_CLOUDFRAC,comm_TOTCLOUDFRAC,& |
---|
| 58 | comm_RAIN,comm_SNOW,comm_ALBEQ,& |
---|
| 59 | comm_FLUXTOP_DN,comm_FLUXABS_SW,& |
---|
| 60 | comm_FLUXTOP_LW,comm_FLUXSURF_SW,& |
---|
| 61 | comm_FLUXSURF_LW,comm_FLXGRD,& |
---|
| 62 | comm_LSCEZ,comm_H2OICE_REFF |
---|
| 63 | #endif |
---|
| 64 | |
---|
[1623] | 65 | #ifdef CPP_XIOS |
---|
[1622] | 66 | use xios_output_mod, only: initialize_xios_output, & |
---|
| 67 | update_xios_timestep, & |
---|
| 68 | send_xios_field |
---|
[1682] | 69 | use wxios, only: wxios_context_init, xios_context_finalize |
---|
[1623] | 70 | #endif |
---|
[1836] | 71 | |
---|
[253] | 72 | implicit none |
---|
| 73 | |
---|
| 74 | |
---|
| 75 | !================================================================== |
---|
| 76 | ! |
---|
| 77 | ! Purpose |
---|
| 78 | ! ------- |
---|
| 79 | ! Central subroutine for all the physics parameterisations in the |
---|
| 80 | ! universal model. Originally adapted from the Mars LMDZ model. |
---|
| 81 | ! |
---|
| 82 | ! The model can be run without or with tracer transport |
---|
| 83 | ! depending on the value of "tracer" in file "callphys.def". |
---|
| 84 | ! |
---|
| 85 | ! |
---|
| 86 | ! It includes: |
---|
| 87 | ! |
---|
[1477] | 88 | ! I. Initialization : |
---|
| 89 | ! I.1 Firstcall initializations. |
---|
| 90 | ! I.2 Initialization for every call to physiq. |
---|
[253] | 91 | ! |
---|
[1477] | 92 | ! II. Compute radiative transfer tendencies (longwave and shortwave) : |
---|
| 93 | ! II.a Option 1 : Call correlated-k radiative transfer scheme. |
---|
| 94 | ! II.b Option 2 : Call Newtonian cooling scheme. |
---|
| 95 | ! II.c Option 3 : Atmosphere has no radiative effect. |
---|
| 96 | ! |
---|
| 97 | ! III. Vertical diffusion (turbulent mixing) : |
---|
| 98 | ! |
---|
| 99 | ! IV. Dry Convective adjustment : |
---|
| 100 | ! |
---|
| 101 | ! V. Condensation and sublimation of gases (currently just CO2). |
---|
| 102 | ! |
---|
| 103 | ! VI. Tracers |
---|
| 104 | ! VI.1. Water and water ice. |
---|
[1801] | 105 | ! VI.2 Photochemistry |
---|
| 106 | ! VI.3. Aerosols and particles. |
---|
| 107 | ! VI.4. Updates (pressure variations, surface budget). |
---|
| 108 | ! VI.5. Slab Ocean. |
---|
| 109 | ! VI.6. Surface Tracer Update. |
---|
[1477] | 110 | ! |
---|
| 111 | ! VII. Surface and sub-surface soil temperature. |
---|
| 112 | ! |
---|
| 113 | ! VIII. Perform diagnostics and write output files. |
---|
| 114 | ! |
---|
| 115 | ! |
---|
[253] | 116 | ! arguments |
---|
| 117 | ! --------- |
---|
| 118 | ! |
---|
[1477] | 119 | ! INPUT |
---|
[253] | 120 | ! ----- |
---|
[1477] | 121 | ! |
---|
[253] | 122 | ! ngrid Size of the horizontal grid. |
---|
| 123 | ! nlayer Number of vertical layers. |
---|
[1477] | 124 | ! nq Number of advected fields. |
---|
| 125 | ! nametrac Name of corresponding advected fields. |
---|
| 126 | ! |
---|
| 127 | ! firstcall True at the first call. |
---|
| 128 | ! lastcall True at the last call. |
---|
| 129 | ! |
---|
| 130 | ! pday Number of days counted from the North. Spring equinoxe. |
---|
| 131 | ! ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT. |
---|
| 132 | ! ptimestep timestep (s). |
---|
| 133 | ! |
---|
| 134 | ! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa). |
---|
| 135 | ! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa). |
---|
| 136 | ! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2.s-2). |
---|
| 137 | ! |
---|
| 138 | ! pu(ngrid,nlayer) u, zonal component of the wind (ms-1). |
---|
| 139 | ! pv(ngrid,nlayer) v, meridional component of the wind (ms-1). |
---|
| 140 | ! |
---|
| 141 | ! pt(ngrid,nlayer) Temperature (K). |
---|
| 142 | ! |
---|
| 143 | ! pq(ngrid,nlayer,nq) Advected fields. |
---|
| 144 | ! |
---|
[1216] | 145 | ! pudyn(ngrid,nlayer) \ |
---|
[253] | 146 | ! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the |
---|
[1477] | 147 | ! ptdyn(ngrid,nlayer) / corresponding variables. |
---|
[253] | 148 | ! pqdyn(ngrid,nlayer,nq) / |
---|
[1312] | 149 | ! flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary |
---|
[253] | 150 | ! |
---|
[1477] | 151 | ! OUTPUT |
---|
[253] | 152 | ! ------ |
---|
| 153 | ! |
---|
[1308] | 154 | ! pdu(ngrid,nlayer) \ |
---|
| 155 | ! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding |
---|
| 156 | ! pdt(ngrid,nlayer) / variables due to physical processes. |
---|
| 157 | ! pdq(ngrid,nlayer) / |
---|
[253] | 158 | ! pdpsrf(ngrid) / |
---|
| 159 | ! |
---|
| 160 | ! |
---|
| 161 | ! Authors |
---|
| 162 | ! ------- |
---|
[1524] | 163 | ! Frederic Hourdin 15/10/93 |
---|
| 164 | ! Francois Forget 1994 |
---|
| 165 | ! Christophe Hourdin 02/1997 |
---|
[253] | 166 | ! Subroutine completely rewritten by F. Forget (01/2000) |
---|
| 167 | ! Water ice clouds: Franck Montmessin (update 06/2003) |
---|
| 168 | ! Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009) |
---|
| 169 | ! New correlated-k radiative scheme: R. Wordsworth (2009) |
---|
| 170 | ! Many specifically Martian subroutines removed: R. Wordsworth (2009) |
---|
| 171 | ! Improved water cycle: R. Wordsworth / B. Charnay (2010) |
---|
| 172 | ! To F90: R. Wordsworth (2010) |
---|
[594] | 173 | ! New turbulent diffusion scheme: J. Leconte (2012) |
---|
[716] | 174 | ! Loops converted to F90 matrix format: J. Leconte (2012) |
---|
[787] | 175 | ! No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012) |
---|
[1477] | 176 | ! Purge of the code : M. Turbet (2015) |
---|
[1801] | 177 | ! Photochemical core developped by F. Lefevre: B. Charnay (2017) |
---|
[253] | 178 | !================================================================== |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | ! 0. Declarations : |
---|
| 182 | ! ------------------ |
---|
| 183 | |
---|
[1669] | 184 | include "netcdf.inc" |
---|
[253] | 185 | |
---|
| 186 | ! Arguments : |
---|
| 187 | ! ----------- |
---|
| 188 | |
---|
[1477] | 189 | ! INPUTS: |
---|
[253] | 190 | ! ------- |
---|
| 191 | |
---|
[1477] | 192 | integer,intent(in) :: ngrid ! Number of atmospheric columns. |
---|
| 193 | integer,intent(in) :: nlayer ! Number of atmospheric layers. |
---|
| 194 | integer,intent(in) :: nq ! Number of tracers. |
---|
| 195 | character*20,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics. |
---|
| 196 | |
---|
| 197 | logical,intent(in) :: firstcall ! Signals first call to physics. |
---|
| 198 | logical,intent(in) :: lastcall ! Signals last call to physics. |
---|
| 199 | |
---|
| 200 | real,intent(in) :: pday ! Number of elapsed sols since reference Ls=0. |
---|
| 201 | real,intent(in) :: ptime ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon). |
---|
| 202 | real,intent(in) :: ptimestep ! Physics timestep (s). |
---|
| 203 | real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa). |
---|
| 204 | real,intent(in) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa). |
---|
| 205 | real,intent(in) :: pphi(ngrid,nlayer) ! Geopotential at mid-layer (m2s-2). |
---|
| 206 | real,intent(in) :: pu(ngrid,nlayer) ! Zonal wind component (m/s). |
---|
| 207 | real,intent(in) :: pv(ngrid,nlayer) ! Meridional wind component (m/s). |
---|
| 208 | real,intent(in) :: pt(ngrid,nlayer) ! Temperature (K). |
---|
| 209 | real,intent(in) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air). |
---|
| 210 | real,intent(in) :: flxw(ngrid,nlayer) ! Vertical mass flux (ks/s) at lower boundary of layer |
---|
[253] | 211 | |
---|
[1477] | 212 | ! OUTPUTS: |
---|
[253] | 213 | ! -------- |
---|
| 214 | |
---|
[1477] | 215 | ! Physical tendencies : |
---|
| 216 | |
---|
| 217 | real,intent(out) :: pdu(ngrid,nlayer) ! Zonal wind tendencies (m/s/s). |
---|
| 218 | real,intent(out) :: pdv(ngrid,nlayer) ! Meridional wind tendencies (m/s/s). |
---|
| 219 | real,intent(out) :: pdt(ngrid,nlayer) ! Temperature tendencies (K/s). |
---|
| 220 | real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s). |
---|
| 221 | real,intent(out) :: pdpsrf(ngrid) ! Surface pressure tendency (Pa/s). |
---|
| 222 | |
---|
[253] | 223 | ! Local saved variables: |
---|
| 224 | ! ---------------------- |
---|
[1622] | 225 | integer,save :: day_ini ! Initial date of the run (sol since Ls=0). |
---|
| 226 | integer,save :: icount ! Counter of calls to physiq during the run. |
---|
| 227 | !$OMP THREADPRIVATE(day_ini,icount) |
---|
| 228 | |
---|
[253] | 229 | ! Local variables : |
---|
| 230 | ! ----------------- |
---|
| 231 | |
---|
[1477] | 232 | ! Aerosol (dust or ice) extinction optical depth at reference wavelength |
---|
[253] | 233 | ! for the "naerkind" optically active aerosols: |
---|
| 234 | |
---|
[1477] | 235 | real aerosol(ngrid,nlayer,naerkind) ! Aerosols. |
---|
| 236 | real zh(ngrid,nlayer) ! Potential temperature (K). |
---|
| 237 | real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!) |
---|
| 238 | |
---|
[1669] | 239 | integer l,ig,ierr,iq,nw,isoil |
---|
[1161] | 240 | |
---|
[1524] | 241 | |
---|
[1477] | 242 | !$OMP zdtlw,zdtsw,sensibFlux) |
---|
[253] | 243 | |
---|
[1477] | 244 | real zls ! Solar longitude (radians). |
---|
| 245 | real zlss ! Sub solar point longitude (radians). |
---|
| 246 | real zday ! Date (time since Ls=0, calculated in sols). |
---|
| 247 | real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers. |
---|
| 248 | real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries. |
---|
[253] | 249 | |
---|
[1477] | 250 | ! TENDENCIES due to various processes : |
---|
[253] | 251 | |
---|
[1477] | 252 | ! For Surface Temperature : (K/s) |
---|
| 253 | real zdtsurf(ngrid) ! Cumulated tendencies. |
---|
| 254 | real zdtsurfmr(ngrid) ! Mass_redistribution routine. |
---|
| 255 | real zdtsurfc(ngrid) ! Condense_co2 routine. |
---|
| 256 | real zdtsdif(ngrid) ! Turbdiff/vdifc routines. |
---|
| 257 | real zdtsurf_hyd(ngrid) ! Hydrol routine. |
---|
| 258 | |
---|
| 259 | ! For Atmospheric Temperatures : (K/s) |
---|
| 260 | real dtlscale(ngrid,nlayer) ! Largescale routine. |
---|
| 261 | real zdtc(ngrid,nlayer) ! Condense_co2 routine. |
---|
[1524] | 262 | real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines. |
---|
[1477] | 263 | real zdtmr(ngrid,nlayer) ! Mass_redistribution routine. |
---|
| 264 | real zdtrain(ngrid,nlayer) ! Rain routine. |
---|
| 265 | real dtmoist(ngrid,nlayer) ! Moistadj routine. |
---|
| 266 | real dt_ekman(ngrid,noceanmx), dt_hdiff(ngrid,noceanmx) ! Slab_ocean routine. |
---|
| 267 | real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine. |
---|
| 268 | |
---|
| 269 | ! For Surface Tracers : (kg/m2/s) |
---|
| 270 | real dqsurf(ngrid,nq) ! Cumulated tendencies. |
---|
[1484] | 271 | real zdqsurfc(ngrid) ! Condense_co2 routine. |
---|
[1477] | 272 | real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines. |
---|
| 273 | real zdqssed(ngrid,nq) ! Callsedim routine. |
---|
| 274 | real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine. |
---|
| 275 | real zdqsrain(ngrid), zdqssnow(ngrid) ! Rain routine. |
---|
| 276 | real dqs_hyd(ngrid,nq) ! Hydrol routine. |
---|
| 277 | |
---|
| 278 | ! For Tracers : (kg/kg_of_air/s) |
---|
| 279 | real zdqc(ngrid,nlayer,nq) ! Condense_co2 routine. |
---|
| 280 | real zdqadj(ngrid,nlayer,nq) ! Convadj routine. |
---|
| 281 | real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines. |
---|
| 282 | real zdqevap(ngrid,nlayer) ! Turbdiff routine. |
---|
| 283 | real zdqsed(ngrid,nlayer,nq) ! Callsedim routine. |
---|
| 284 | real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine. |
---|
| 285 | real zdqrain(ngrid,nlayer,nq) ! Rain routine. |
---|
| 286 | real dqmoist(ngrid,nlayer,nq) ! Moistadj routine. |
---|
| 287 | real dqvaplscale(ngrid,nlayer) ! Largescale routine. |
---|
| 288 | real dqcldlscale(ngrid,nlayer) ! Largescale routine. |
---|
[1801] | 289 | REAL zdqchim(ngrid,nlayer,nq) ! Calchim_asis routine |
---|
| 290 | REAL zdqschim(ngrid,nq) ! Calchim_asis routine |
---|
| 291 | |
---|
| 292 | REAL array_zero1(ngrid) |
---|
| 293 | REAL array_zero2(ngrid,nlayer) |
---|
[1477] | 294 | |
---|
| 295 | ! For Winds : (m/s/s) |
---|
| 296 | real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine. |
---|
| 297 | real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer) ! Mass_redistribution routine. |
---|
| 298 | real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines. |
---|
| 299 | real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines. |
---|
| 300 | real zdhadj(ngrid,nlayer) ! Convadj routine. |
---|
| 301 | |
---|
| 302 | ! For Pressure and Mass : |
---|
| 303 | real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s). |
---|
| 304 | real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s). |
---|
| 305 | real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s). |
---|
[253] | 306 | |
---|
[1477] | 307 | |
---|
| 308 | |
---|
| 309 | ! Local variables for LOCAL CALCULATIONS: |
---|
| 310 | ! --------------------------------------- |
---|
[787] | 311 | real zflubid(ngrid) |
---|
[1308] | 312 | real zplanck(ngrid),zpopsk(ngrid,nlayer) |
---|
[253] | 313 | real ztim1,ztim2,ztim3, z1,z2 |
---|
| 314 | real ztime_fin |
---|
[1308] | 315 | real zdh(ngrid,nlayer) |
---|
[1194] | 316 | real gmplanet |
---|
[1297] | 317 | real taux(ngrid),tauy(ngrid) |
---|
[1194] | 318 | |
---|
[253] | 319 | |
---|
[1477] | 320 | ! local variables for DIAGNOSTICS : (diagfi & stat) |
---|
| 321 | ! ------------------------------------------------- |
---|
| 322 | real ps(ngrid) ! Surface Pressure. |
---|
| 323 | real zt(ngrid,nlayer) ! Atmospheric Temperature. |
---|
| 324 | real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds. |
---|
| 325 | real zq(ngrid,nlayer,nq) ! Atmospheric Tracers. |
---|
| 326 | real zdtadj(ngrid,nlayer) ! Convadj Diagnostic. |
---|
| 327 | real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s). |
---|
[1637] | 328 | real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2). |
---|
[253] | 329 | |
---|
[1477] | 330 | real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T). |
---|
| 331 | real vmr(ngrid,nlayer) ! volume mixing ratio |
---|
[253] | 332 | real time_phys |
---|
[597] | 333 | |
---|
[1477] | 334 | real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic. |
---|
| 335 | |
---|
| 336 | real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2). |
---|
[253] | 337 | |
---|
| 338 | ! included by RW for H2O Manabe scheme |
---|
[1477] | 339 | real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj). |
---|
| 340 | real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale). |
---|
[253] | 341 | |
---|
| 342 | |
---|
[594] | 343 | ! to test energy conservation (RW+JL) |
---|
[1308] | 344 | real mass(ngrid,nlayer),massarea(ngrid,nlayer) |
---|
[651] | 345 | real dEtot, dEtots, AtmToSurf_TurbFlux |
---|
[959] | 346 | real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW |
---|
[1315] | 347 | !$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW) |
---|
[1308] | 348 | real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer) |
---|
[787] | 349 | real dEdiffs(ngrid),dEdiff(ngrid) |
---|
[1308] | 350 | real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayer), lscaledEz(ngrid,nlayer) |
---|
[1477] | 351 | |
---|
[594] | 352 | !JL12 conservation test for mean flow kinetic energy has been disabled temporarily |
---|
[1477] | 353 | |
---|
| 354 | real dtmoist_max,dtmoist_min |
---|
[1295] | 355 | real dItot, dItot_tmp, dVtot, dVtot_tmp |
---|
[253] | 356 | |
---|
| 357 | |
---|
[1477] | 358 | real h2otot ! Total Amount of water. For diagnostic. |
---|
| 359 | real icesrf,liqsrf,icecol,vapcol ! Total Amounts of water (ice,liq,vap). For diagnostic. |
---|
[1295] | 360 | real dWtot, dWtot_tmp, dWtots, dWtots_tmp |
---|
[1477] | 361 | logical,save :: watertest |
---|
[1315] | 362 | !$OMP THREADPRIVATE(watertest) |
---|
[253] | 363 | |
---|
[1477] | 364 | real qsat(ngrid,nlayer) ! Water Vapor Volume Mixing Ratio at saturation (kg/kg_of_air). |
---|
| 365 | real RH(ngrid,nlayer) ! Relative humidity. |
---|
| 366 | real H2Omaxcol(ngrid) ! Maximum possible H2O column amount (at 100% saturation) (kg/m2). |
---|
| 367 | real psat_tmp |
---|
| 368 | |
---|
| 369 | logical clearsky ! For double radiative transfer call. By BC |
---|
| 370 | |
---|
[1482] | 371 | ! For Clear Sky Case. |
---|
| 372 | real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid) ! For SW/LW flux. |
---|
| 373 | real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux. |
---|
| 374 | real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation. |
---|
| 375 | real tau_col1(ngrid) ! For aerosol optical depth diagnostic. |
---|
| 376 | real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV) ! For Outgoing Radiation diagnostics. |
---|
[253] | 377 | real tf, ntf |
---|
| 378 | |
---|
[1477] | 379 | real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW |
---|
[253] | 380 | |
---|
[1477] | 381 | real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW |
---|
[253] | 382 | |
---|
[787] | 383 | real reffcol(ngrid,naerkind) |
---|
[253] | 384 | |
---|
[1477] | 385 | ! Sourceevol for 'accelerated ice evolution'. By RW |
---|
[305] | 386 | real delta_ice,ice_tot |
---|
[253] | 387 | integer num_run |
---|
[728] | 388 | logical,save :: ice_update |
---|
[996] | 389 | |
---|
[1297] | 390 | |
---|
| 391 | real :: tsurf2(ngrid) |
---|
| 392 | real :: flux_o(ngrid),flux_g(ngrid),fluxgrdocean(ngrid) |
---|
| 393 | real :: flux_sens_lat(ngrid) |
---|
| 394 | real :: qsurfint(ngrid,nq) |
---|
| 395 | |
---|
| 396 | |
---|
[1477] | 397 | !================================================================================================== |
---|
[253] | 398 | |
---|
| 399 | ! ----------------- |
---|
[1477] | 400 | ! I. INITIALISATION |
---|
| 401 | ! ----------------- |
---|
[253] | 402 | |
---|
[1477] | 403 | ! -------------------------------- |
---|
| 404 | ! I.1 First Call Initialisation. |
---|
| 405 | ! -------------------------------- |
---|
[253] | 406 | if (firstcall) then |
---|
[1477] | 407 | ! Allocate saved arrays. |
---|
[1836] | 408 | call phys_stat_var_mod |
---|
[858] | 409 | |
---|
[1477] | 410 | ! Variables set to 0 |
---|
[253] | 411 | ! ~~~~~~~~~~~~~~~~~~ |
---|
| 412 | dtrad(:,:) = 0.0 |
---|
| 413 | fluxrad(:) = 0.0 |
---|
| 414 | tau_col(:) = 0.0 |
---|
| 415 | zdtsw(:,:) = 0.0 |
---|
| 416 | zdtlw(:,:) = 0.0 |
---|
[726] | 417 | |
---|
| 418 | |
---|
[1477] | 419 | ! Initialize aerosol indexes. |
---|
| 420 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 421 | call iniaerosol() |
---|
| 422 | |
---|
[253] | 423 | |
---|
[1477] | 424 | ! Initialize tracer names, indexes and properties. |
---|
| 425 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 426 | IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) |
---|
[253] | 427 | if (tracer) then |
---|
[787] | 428 | call initracer(ngrid,nq,nametrac) |
---|
[1801] | 429 | if(photochem) then |
---|
| 430 | call ini_conc_mod(ngrid,nlayer) |
---|
| 431 | endif |
---|
[1477] | 432 | endif |
---|
[253] | 433 | |
---|
[1682] | 434 | #ifdef CPP_XIOS |
---|
| 435 | ! Initialize XIOS context |
---|
| 436 | write(*,*) "physiq: call wxios_context_init" |
---|
| 437 | CALL wxios_context_init |
---|
| 438 | #endif |
---|
[726] | 439 | |
---|
[1477] | 440 | ! Read 'startfi.nc' file. |
---|
[253] | 441 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1836] | 442 | #ifndef MESOSCALE |
---|
[1669] | 443 | call phyetat0(startphy_file, & |
---|
| 444 | ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, & |
---|
[1477] | 445 | day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, & |
---|
| 446 | cloudfrac,totcloudfrac,hice, & |
---|
| 447 | rnat,pctsrf_sic,tslab, tsea_ice,sea_ice) |
---|
[1836] | 448 | #else |
---|
| 449 | emis(:)=0.0 |
---|
| 450 | q2(:,:)=0.0 |
---|
| 451 | qsurf(:,:)=0.0 |
---|
| 452 | day_ini = pday |
---|
| 453 | #endif |
---|
| 454 | |
---|
| 455 | #ifndef MESOSCALE |
---|
[1669] | 456 | if (.not.startphy_file) then |
---|
| 457 | ! additionnal "academic" initialization of physics |
---|
| 458 | if (is_master) write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!" |
---|
| 459 | tsurf(:)=pt(:,1) |
---|
| 460 | if (is_master) write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!" |
---|
| 461 | do isoil=1,nsoilmx |
---|
| 462 | tsoil(1:ngrid,isoil)=tsurf(1:ngrid) |
---|
| 463 | enddo |
---|
| 464 | if (is_master) write(*,*) "Physiq: initializing day_ini to pdat !" |
---|
| 465 | day_ini=pday |
---|
| 466 | endif |
---|
[1836] | 467 | #endif |
---|
[253] | 468 | |
---|
| 469 | if (pday.ne.day_ini) then |
---|
| 470 | write(*,*) "ERROR in physiq.F90:" |
---|
| 471 | write(*,*) "bad synchronization between physics and dynamics" |
---|
| 472 | write(*,*) "dynamics day: ",pday |
---|
| 473 | write(*,*) "physics day: ",day_ini |
---|
| 474 | stop |
---|
| 475 | endif |
---|
| 476 | |
---|
| 477 | write (*,*) 'In physiq day_ini =', day_ini |
---|
| 478 | |
---|
[1482] | 479 | ! Initialize albedo calculation. |
---|
| 480 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 481 | albedo(:,:)=0.0 |
---|
[1524] | 482 | albedo_bareground(:)=0.0 |
---|
| 483 | albedo_snow_SPECTV(:)=0.0 |
---|
| 484 | albedo_co2_ice_SPECTV(:)=0.0 |
---|
[1482] | 485 | call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) |
---|
| 486 | |
---|
| 487 | ! Initialize orbital calculation. |
---|
| 488 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[253] | 489 | call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) |
---|
| 490 | |
---|
| 491 | |
---|
| 492 | if(tlocked)then |
---|
| 493 | print*,'Planet is tidally locked at resonance n=',nres |
---|
| 494 | print*,'Make sure you have the right rotation rate!!!' |
---|
| 495 | endif |
---|
| 496 | |
---|
[1477] | 497 | ! Initialize oceanic variables. |
---|
| 498 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1297] | 499 | |
---|
| 500 | if (ok_slab_ocean)then |
---|
| 501 | |
---|
| 502 | call ocean_slab_init(ngrid,ptimestep, tslab, & |
---|
[1477] | 503 | sea_ice, pctsrf_sic) |
---|
[1297] | 504 | |
---|
[1529] | 505 | call ini_surf_heat_transp_mod() |
---|
| 506 | |
---|
[1477] | 507 | knindex(:) = 0 |
---|
[1297] | 508 | |
---|
[1477] | 509 | do ig=1,ngrid |
---|
[1297] | 510 | zmasq(ig)=1 |
---|
| 511 | knindex(ig) = ig |
---|
| 512 | if (nint(rnat(ig)).eq.0) then |
---|
| 513 | zmasq(ig)=0 |
---|
| 514 | endif |
---|
[1477] | 515 | enddo |
---|
[1297] | 516 | |
---|
[1308] | 517 | CALL init_masquv(ngrid,zmasq) |
---|
[1297] | 518 | |
---|
[1477] | 519 | endif ! end of 'ok_slab_ocean'. |
---|
[1297] | 520 | |
---|
| 521 | |
---|
[1477] | 522 | ! Initialize soil. |
---|
| 523 | ! ~~~~~~~~~~~~~~~~ |
---|
[253] | 524 | if (callsoil) then |
---|
[1477] | 525 | |
---|
[787] | 526 | call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, & |
---|
[1477] | 527 | ptimestep,tsurf,tsoil,capcal,fluxgrd) |
---|
[1297] | 528 | |
---|
| 529 | if (ok_slab_ocean) then |
---|
| 530 | do ig=1,ngrid |
---|
| 531 | if (nint(rnat(ig)).eq.2) then |
---|
[1477] | 532 | capcal(ig)=capcalocean |
---|
| 533 | if (pctsrf_sic(ig).gt.0.5) then |
---|
| 534 | capcal(ig)=capcalseaice |
---|
| 535 | if (qsurf(ig,igcm_h2o_ice).gt.0.) then |
---|
| 536 | capcal(ig)=capcalsno |
---|
| 537 | endif |
---|
| 538 | endif |
---|
[1297] | 539 | endif |
---|
| 540 | enddo |
---|
[1477] | 541 | endif ! end of 'ok_slab_ocean'. |
---|
[1297] | 542 | |
---|
[1477] | 543 | else ! else of 'callsoil'. |
---|
| 544 | |
---|
[253] | 545 | print*,'WARNING! Thermal conduction in the soil turned off' |
---|
[918] | 546 | capcal(:)=1.e6 |
---|
[952] | 547 | fluxgrd(:)=intheat |
---|
| 548 | print*,'Flux from ground = ',intheat,' W m^-2' |
---|
[1477] | 549 | |
---|
| 550 | endif ! end of 'callsoil'. |
---|
| 551 | |
---|
[253] | 552 | icount=1 |
---|
| 553 | |
---|
[1477] | 554 | ! Decide whether to update ice at end of run. |
---|
| 555 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[253] | 556 | ice_update=.false. |
---|
| 557 | if(sourceevol)then |
---|
[1315] | 558 | !$OMP MASTER |
---|
[955] | 559 | open(128,file='num_run',form='formatted', & |
---|
| 560 | status="old",iostat=ierr) |
---|
| 561 | if (ierr.ne.0) then |
---|
[1477] | 562 | write(*,*) "physiq: Error! No num_run file!" |
---|
| 563 | write(*,*) " (which is needed for sourceevol option)" |
---|
| 564 | stop |
---|
[955] | 565 | endif |
---|
[253] | 566 | read(128,*) num_run |
---|
| 567 | close(128) |
---|
[1315] | 568 | !$OMP END MASTER |
---|
| 569 | !$OMP BARRIER |
---|
[253] | 570 | |
---|
[365] | 571 | if(num_run.ne.0.and.mod(num_run,2).eq.0)then |
---|
[253] | 572 | print*,'Updating ice at end of this year!' |
---|
| 573 | ice_update=.true. |
---|
| 574 | ice_min(:)=1.e4 |
---|
[1477] | 575 | endif |
---|
| 576 | |
---|
| 577 | endif ! end of 'sourceevol'. |
---|
[253] | 578 | |
---|
[1477] | 579 | |
---|
| 580 | ! Here is defined the type of the surface : Continent or Ocean. |
---|
| 581 | ! BC2014 : This is now already done in newstart. |
---|
| 582 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1297] | 583 | if (.not.ok_slab_ocean) then |
---|
[1477] | 584 | |
---|
[1297] | 585 | rnat(:)=1. |
---|
| 586 | do ig=1,ngrid |
---|
[1477] | 587 | if(inertiedat(ig,1).gt.1.E4)then |
---|
| 588 | rnat(ig)=0 |
---|
| 589 | endif |
---|
[1297] | 590 | enddo |
---|
[253] | 591 | |
---|
[1297] | 592 | print*,'WARNING! Surface type currently decided by surface inertia' |
---|
| 593 | print*,'This should be improved e.g. in newstart.F' |
---|
[1477] | 594 | |
---|
| 595 | endif ! end of 'ok_slab_ocean'. |
---|
[253] | 596 | |
---|
[1477] | 597 | |
---|
| 598 | ! Initialize surface history variable. |
---|
[253] | 599 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[728] | 600 | qsurf_hist(:,:)=qsurf(:,:) |
---|
[253] | 601 | |
---|
[1637] | 602 | ! Initialize variable for dynamical heating and zonal wind tendency diagnostic |
---|
| 603 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[253] | 604 | ztprevious(:,:)=pt(:,:) |
---|
[1637] | 605 | zuprevious(:,:)=pu(:,:) |
---|
[253] | 606 | |
---|
| 607 | ! Set temperature just above condensation temperature (for Early Mars) |
---|
| 608 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 609 | if(nearco2cond) then |
---|
| 610 | write(*,*)' WARNING! Starting at Tcond+1K' |
---|
| 611 | do l=1, nlayer |
---|
| 612 | do ig=1,ngrid |
---|
| 613 | pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 & |
---|
| 614 | -pt(ig,l)) / ptimestep |
---|
| 615 | enddo |
---|
| 616 | enddo |
---|
| 617 | endif |
---|
| 618 | |
---|
[1477] | 619 | if(meanOLR)then |
---|
| 620 | call system('rm -f rad_bal.out') ! to record global radiative balance. |
---|
| 621 | call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures. |
---|
| 622 | call system('rm -f h2o_bal.out') ! to record global hydrological balance. |
---|
[253] | 623 | endif |
---|
| 624 | |
---|
| 625 | |
---|
[1477] | 626 | watertest=.false. |
---|
| 627 | if(water)then ! Initialize variables for water cycle |
---|
| 628 | |
---|
[365] | 629 | if(enertest)then |
---|
| 630 | watertest = .true. |
---|
| 631 | endif |
---|
| 632 | |
---|
[728] | 633 | if(ice_update)then |
---|
| 634 | ice_initial(:)=qsurf(:,igcm_h2o_ice) |
---|
| 635 | endif |
---|
[253] | 636 | |
---|
| 637 | endif |
---|
[1477] | 638 | |
---|
[253] | 639 | call su_watercycle ! even if we don't have a water cycle, we might |
---|
| 640 | ! need epsi for the wvp definitions in callcorrk.F |
---|
[1836] | 641 | #ifndef MESOSCALE |
---|
[1477] | 642 | if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d. |
---|
[1542] | 643 | call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & |
---|
| 644 | ptimestep,pday+nday,time_phys,cell_area, & |
---|
[1482] | 645 | albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) |
---|
[1216] | 646 | endif |
---|
[1836] | 647 | #endif |
---|
[1216] | 648 | |
---|
[1622] | 649 | ! XIOS outputs |
---|
| 650 | #ifdef CPP_XIOS |
---|
| 651 | |
---|
| 652 | write(*,*) "physiq: call initialize_xios_output" |
---|
| 653 | call initialize_xios_output(pday,ptime,ptimestep,daysec, & |
---|
| 654 | presnivs,pseudoalt) |
---|
| 655 | #endif |
---|
[1682] | 656 | write(*,*) "physiq: end of firstcall" |
---|
[1477] | 657 | endif ! end of 'firstcall' |
---|
[253] | 658 | |
---|
[1477] | 659 | ! ------------------------------------------------------ |
---|
| 660 | ! I.2 Initializations done at every physical timestep: |
---|
| 661 | ! ------------------------------------------------------ |
---|
| 662 | |
---|
[1622] | 663 | #ifdef CPP_XIOS |
---|
| 664 | ! update XIOS time/calendar |
---|
| 665 | call update_xios_timestep |
---|
| 666 | #endif |
---|
| 667 | |
---|
[1477] | 668 | ! Initialize various variables |
---|
| 669 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1297] | 670 | |
---|
[253] | 671 | if ( .not.nearco2cond ) then |
---|
[1308] | 672 | pdt(1:ngrid,1:nlayer) = 0.0 |
---|
[1477] | 673 | endif |
---|
| 674 | zdtsurf(1:ngrid) = 0.0 |
---|
[1308] | 675 | pdq(1:ngrid,1:nlayer,1:nq) = 0.0 |
---|
[1477] | 676 | dqsurf(1:ngrid,1:nq)= 0.0 |
---|
| 677 | pdu(1:ngrid,1:nlayer) = 0.0 |
---|
| 678 | pdv(1:ngrid,1:nlayer) = 0.0 |
---|
[787] | 679 | pdpsrf(1:ngrid) = 0.0 |
---|
[1477] | 680 | zflubid(1:ngrid) = 0.0 |
---|
[1297] | 681 | flux_sens_lat(1:ngrid) = 0.0 |
---|
| 682 | taux(1:ngrid) = 0.0 |
---|
| 683 | tauy(1:ngrid) = 0.0 |
---|
[253] | 684 | |
---|
[1477] | 685 | zday=pday+ptime ! Compute time, in sols (and fraction thereof). |
---|
[1297] | 686 | |
---|
[1477] | 687 | ! Compute Stellar Longitude (Ls), and orbital parameters. |
---|
| 688 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[253] | 689 | if (season) then |
---|
| 690 | call stellarlong(zday,zls) |
---|
| 691 | else |
---|
| 692 | call stellarlong(float(day_ini),zls) |
---|
| 693 | end if |
---|
| 694 | |
---|
[1329] | 695 | call orbite(zls,dist_star,declin,right_ascen) |
---|
[1477] | 696 | |
---|
[1329] | 697 | if (tlocked) then |
---|
| 698 | zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi) |
---|
| 699 | elseif (diurnal) then |
---|
[1524] | 700 | zlss=-2.*pi*(zday-.5) |
---|
[1329] | 701 | else if(diurnal .eqv. .false.) then |
---|
| 702 | zlss=9999. |
---|
| 703 | endif |
---|
[1194] | 704 | |
---|
| 705 | |
---|
[1477] | 706 | ! Compute variations of g with latitude (oblate case). |
---|
| 707 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 708 | if (oblate .eqv. .false.) then |
---|
| 709 | glat(:) = g |
---|
| 710 | else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then |
---|
| 711 | print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)' |
---|
| 712 | call abort |
---|
| 713 | else |
---|
| 714 | gmplanet = MassPlanet*grav*1e24 |
---|
| 715 | do ig=1,ngrid |
---|
[1542] | 716 | glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - latitude(ig)))) |
---|
[1477] | 717 | end do |
---|
| 718 | endif |
---|
[1297] | 719 | |
---|
[1477] | 720 | ! Compute geopotential between layers. |
---|
| 721 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1308] | 722 | zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer) |
---|
| 723 | do l=1,nlayer |
---|
[1477] | 724 | zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid) |
---|
[1194] | 725 | enddo |
---|
[728] | 726 | |
---|
[787] | 727 | zzlev(1:ngrid,1)=0. |
---|
[1477] | 728 | zzlev(1:ngrid,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km... |
---|
[728] | 729 | |
---|
[253] | 730 | do l=2,nlayer |
---|
| 731 | do ig=1,ngrid |
---|
| 732 | z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) |
---|
| 733 | z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) |
---|
| 734 | zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) |
---|
| 735 | enddo |
---|
[1477] | 736 | enddo |
---|
[253] | 737 | |
---|
[1477] | 738 | ! Compute potential temperature |
---|
| 739 | ! Note : Potential temperature calculation may not be the same in physiq and dynamic... |
---|
| 740 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[597] | 741 | do l=1,nlayer |
---|
[787] | 742 | do ig=1,ngrid |
---|
[253] | 743 | zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp |
---|
[597] | 744 | zh(ig,l)=pt(ig,l)/zpopsk(ig,l) |
---|
[1194] | 745 | mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig) |
---|
[1542] | 746 | massarea(ig,l)=mass(ig,l)*cell_area(ig) |
---|
[253] | 747 | enddo |
---|
| 748 | enddo |
---|
| 749 | |
---|
[1312] | 750 | ! Compute vertical velocity (m/s) from vertical mass flux |
---|
[1346] | 751 | ! w = F / (rho*area) and rho = P/(r*T) |
---|
[1477] | 752 | ! But first linearly interpolate mass flux to mid-layers |
---|
| 753 | do l=1,nlayer-1 |
---|
| 754 | pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1)) |
---|
| 755 | enddo |
---|
| 756 | pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0 |
---|
| 757 | do l=1,nlayer |
---|
| 758 | pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) / & |
---|
[1542] | 759 | (pplay(1:ngrid,l)*cell_area(1:ngrid)) |
---|
[1477] | 760 | enddo |
---|
[1194] | 761 | |
---|
[1801] | 762 | ! ---------------------------------------------------------------- |
---|
| 763 | ! Compute mean mass, cp, and R |
---|
| 764 | ! -------------------------------- |
---|
| 765 | |
---|
| 766 | if(photochem) then |
---|
| 767 | call concentrations(ngrid,nlayer,nq,pplay,pt,pdt,pq,pdq,ptimestep) |
---|
| 768 | endif |
---|
| 769 | |
---|
[1477] | 770 | !--------------------------------- |
---|
| 771 | ! II. Compute radiative tendencies |
---|
| 772 | !--------------------------------- |
---|
[253] | 773 | |
---|
| 774 | if (callrad) then |
---|
[526] | 775 | if( mod(icount-1,iradia).eq.0.or.lastcall) then |
---|
[253] | 776 | |
---|
[1477] | 777 | ! Compute local stellar zenith angles |
---|
| 778 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 779 | if (tlocked) then |
---|
| 780 | ! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb |
---|
| 781 | ztim1=SIN(declin) |
---|
| 782 | ztim2=COS(declin)*COS(zlss) |
---|
| 783 | ztim3=COS(declin)*SIN(zlss) |
---|
[253] | 784 | |
---|
[1477] | 785 | call stelang(ngrid,sinlon,coslon,sinlat,coslat, & |
---|
| 786 | ztim1,ztim2,ztim3,mu0,fract, flatten) |
---|
[253] | 787 | |
---|
[1477] | 788 | elseif (diurnal) then |
---|
| 789 | ztim1=SIN(declin) |
---|
| 790 | ztim2=COS(declin)*COS(2.*pi*(zday-.5)) |
---|
| 791 | ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) |
---|
[253] | 792 | |
---|
[1477] | 793 | call stelang(ngrid,sinlon,coslon,sinlat,coslat, & |
---|
| 794 | ztim1,ztim2,ztim3,mu0,fract, flatten) |
---|
| 795 | else if(diurnal .eqv. .false.) then |
---|
[253] | 796 | |
---|
[1542] | 797 | call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten) |
---|
[1161] | 798 | ! WARNING: this function appears not to work in 1D |
---|
[253] | 799 | |
---|
[1477] | 800 | endif |
---|
[1161] | 801 | |
---|
[1477] | 802 | ! Eclipse incoming sunlight (e.g. Saturn ring shadowing). |
---|
[1429] | 803 | if(rings_shadow) then |
---|
| 804 | call call_rings(ngrid, ptime, pday, diurnal) |
---|
| 805 | endif |
---|
[1133] | 806 | |
---|
[1329] | 807 | |
---|
[1477] | 808 | if (corrk) then |
---|
| 809 | |
---|
| 810 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 811 | ! II.a Call correlated-k radiative transfer scheme |
---|
| 812 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 813 | if(kastprof)then |
---|
| 814 | print*,'kastprof should not = true here' |
---|
| 815 | call abort |
---|
| 816 | endif |
---|
[1524] | 817 | if(water) then |
---|
[1308] | 818 | muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap)) |
---|
[1524] | 819 | muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) |
---|
| 820 | ! take into account water effect on mean molecular weight |
---|
| 821 | else |
---|
[1308] | 822 | muvar(1:ngrid,1:nlayer+1)=mugaz |
---|
[1524] | 823 | endif |
---|
[538] | 824 | |
---|
[1297] | 825 | |
---|
[1477] | 826 | if(ok_slab_ocean) then |
---|
| 827 | tsurf2(:)=tsurf(:) |
---|
| 828 | do ig=1,ngrid |
---|
| 829 | if (nint(rnat(ig))==0) then |
---|
| 830 | tsurf(ig)=((1.-pctsrf_sic(ig))*tslab(ig,1)**4+pctsrf_sic(ig)*tsea_ice(ig)**4)**0.25 |
---|
| 831 | endif |
---|
| 832 | enddo |
---|
| 833 | endif !(ok_slab_ocean) |
---|
[1297] | 834 | |
---|
[1477] | 835 | ! standard callcorrk |
---|
| 836 | clearsky=.false. |
---|
[1482] | 837 | call callcorrk(ngrid,nlayer,pq,nq,qsurf, & |
---|
| 838 | albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & |
---|
| 839 | tsurf,fract,dist_star,aerosol,muvar, & |
---|
| 840 | zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, & |
---|
| 841 | fluxsurfabs_sw,fluxtop_lw, & |
---|
| 842 | fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, & |
---|
| 843 | tau_col,cloudfrac,totcloudfrac, & |
---|
| 844 | clearsky,firstcall,lastcall) |
---|
[1297] | 845 | |
---|
[1482] | 846 | ! Option to call scheme once more for clear regions |
---|
[1477] | 847 | if(CLFvarying)then |
---|
[253] | 848 | |
---|
[1477] | 849 | ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ... |
---|
| 850 | clearsky=.true. |
---|
[1482] | 851 | call callcorrk(ngrid,nlayer,pq,nq,qsurf, & |
---|
| 852 | albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt, & |
---|
| 853 | tsurf,fract,dist_star,aerosol,muvar, & |
---|
| 854 | zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1, & |
---|
| 855 | fluxsurfabs_sw1,fluxtop_lw1, & |
---|
| 856 | fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, & |
---|
| 857 | tau_col1,cloudfrac,totcloudfrac, & |
---|
[1477] | 858 | clearsky,firstcall,lastcall) |
---|
| 859 | clearsky = .false. ! just in case. |
---|
[253] | 860 | |
---|
[1477] | 861 | ! Sum the fluxes and heating rates from cloudy/clear cases |
---|
| 862 | do ig=1,ngrid |
---|
| 863 | tf=totcloudfrac(ig) |
---|
[1482] | 864 | ntf=1.-tf |
---|
| 865 | fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig) |
---|
| 866 | fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig) |
---|
| 867 | albedo_equivalent(ig) = ntf*albedo_equivalent1(ig) + tf*albedo_equivalent(ig) |
---|
| 868 | fluxsurfabs_sw(ig) = ntf*fluxsurfabs_sw1(ig) + tf*fluxsurfabs_sw(ig) |
---|
| 869 | fluxtop_lw(ig) = ntf*fluxtop_lw1(ig) + tf*fluxtop_lw(ig) |
---|
| 870 | fluxabs_sw(ig) = ntf*fluxabs_sw1(ig) + tf*fluxabs_sw(ig) |
---|
| 871 | tau_col(ig) = ntf*tau_col1(ig) + tf*tau_col(ig) |
---|
[253] | 872 | |
---|
[1477] | 873 | zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer) |
---|
| 874 | zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer) |
---|
[253] | 875 | |
---|
[1524] | 876 | OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV) |
---|
| 877 | OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI) |
---|
[1482] | 878 | enddo |
---|
[253] | 879 | |
---|
[1477] | 880 | endif ! end of CLFvarying. |
---|
[253] | 881 | |
---|
[1477] | 882 | if(ok_slab_ocean) then |
---|
| 883 | tsurf(:)=tsurf2(:) |
---|
| 884 | endif |
---|
[1297] | 885 | |
---|
| 886 | |
---|
[1482] | 887 | ! Radiative flux from the sky absorbed by the surface (W.m-2). |
---|
[1477] | 888 | GSR=0.0 |
---|
[1482] | 889 | fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid) |
---|
[253] | 890 | |
---|
[1477] | 891 | !if(noradsurf)then ! no lower surface; SW flux just disappears |
---|
[1542] | 892 | ! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea |
---|
[1477] | 893 | ! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid) |
---|
| 894 | ! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2' |
---|
| 895 | !endif |
---|
[253] | 896 | |
---|
[1477] | 897 | ! Net atmospheric radiative heating rate (K.s-1) |
---|
| 898 | dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer) |
---|
[1498] | 899 | |
---|
| 900 | ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that ! |
---|
| 901 | if (firstcall .and. albedo_spectral_mode) then |
---|
| 902 | call spectral_albedo_calc(albedo_snow_SPECTV) |
---|
| 903 | endif |
---|
[253] | 904 | |
---|
[1477] | 905 | elseif(newtonian)then |
---|
[1482] | 906 | |
---|
| 907 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 908 | ! II.b Call Newtonian cooling scheme |
---|
| 909 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1477] | 910 | call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall) |
---|
[253] | 911 | |
---|
[1477] | 912 | zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep |
---|
| 913 | ! e.g. surface becomes proxy for 1st atmospheric layer ? |
---|
[253] | 914 | |
---|
[1477] | 915 | else |
---|
| 916 | |
---|
| 917 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 918 | ! II.c Atmosphere has no radiative effect |
---|
| 919 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 920 | fluxtop_dn(1:ngrid) = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2 |
---|
| 921 | if(ngrid.eq.1)then ! / by 4 globally in 1D case... |
---|
| 922 | fluxtop_dn(1) = fract(1)*Fat1AU/dist_star**2/2.0 |
---|
| 923 | endif |
---|
| 924 | fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid) |
---|
[1482] | 925 | print*,'------------WARNING---WARNING------------' ! by MT2015. |
---|
| 926 | print*,'You are in corrk=false mode, ' |
---|
[1498] | 927 | print*,'and the surface albedo is taken equal to the first visible spectral value' |
---|
| 928 | |
---|
| 929 | fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1)) |
---|
| 930 | fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid) |
---|
[1477] | 931 | fluxtop_lw(1:ngrid) = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4 |
---|
[253] | 932 | |
---|
[1477] | 933 | dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating |
---|
[253] | 934 | |
---|
[1477] | 935 | endif ! end of corrk |
---|
[253] | 936 | |
---|
[1477] | 937 | endif ! of if(mod(icount-1,iradia).eq.0) |
---|
[787] | 938 | |
---|
[253] | 939 | |
---|
[1477] | 940 | ! Transformation of the radiative tendencies |
---|
| 941 | ! ------------------------------------------ |
---|
| 942 | zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid) |
---|
| 943 | zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid) |
---|
| 944 | fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid) |
---|
| 945 | pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer) |
---|
| 946 | |
---|
| 947 | ! Test of energy conservation |
---|
| 948 | !---------------------------- |
---|
[253] | 949 | if(enertest)then |
---|
[1524] | 950 | call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW) |
---|
| 951 | call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW) |
---|
[1542] | 952 | !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk |
---|
| 953 | call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk |
---|
| 954 | call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*cell_area(:)/totarea_planet,dEtotsLW) |
---|
[1524] | 955 | dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:) |
---|
| 956 | dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:) |
---|
| 957 | if (is_master) then |
---|
[1477] | 958 | print*,'---------------------------------------------------------------' |
---|
| 959 | print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2' |
---|
| 960 | print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2' |
---|
| 961 | print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2' |
---|
| 962 | print*,'In corrk SW surface heating =',dEtotsSW,' W m-2' |
---|
| 963 | print*,'In corrk LW surface heating =',dEtotsLW,' W m-2' |
---|
| 964 | print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2' |
---|
[1524] | 965 | endif |
---|
[1477] | 966 | endif ! end of 'enertest' |
---|
[253] | 967 | |
---|
| 968 | endif ! of if (callrad) |
---|
| 969 | |
---|
| 970 | |
---|
[1477] | 971 | |
---|
| 972 | ! -------------------------------------------- |
---|
| 973 | ! III. Vertical diffusion (turbulent mixing) : |
---|
| 974 | ! -------------------------------------------- |
---|
| 975 | |
---|
[253] | 976 | if (calldifv) then |
---|
[526] | 977 | |
---|
[787] | 978 | zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid) |
---|
[253] | 979 | |
---|
[1477] | 980 | ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff. |
---|
[1524] | 981 | if (UseTurbDiff) then |
---|
| 982 | |
---|
[1477] | 983 | call turbdiff(ngrid,nlayer,nq,rnat, & |
---|
| 984 | ptimestep,capcal,lwrite, & |
---|
| 985 | pplay,pplev,zzlay,zzlev,z0, & |
---|
| 986 | pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, & |
---|
| 987 | pdt,pdq,zflubid, & |
---|
| 988 | zdudif,zdvdif,zdtdif,zdtsdif, & |
---|
[1524] | 989 | sensibFlux,q2,zdqdif,zdqevap,zdqsdif, & |
---|
[1477] | 990 | taux,tauy,lastcall) |
---|
[594] | 991 | |
---|
[1524] | 992 | else |
---|
| 993 | |
---|
[1477] | 994 | zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer) |
---|
[594] | 995 | |
---|
[1477] | 996 | call vdifc(ngrid,nlayer,nq,rnat,zpopsk, & |
---|
| 997 | ptimestep,capcal,lwrite, & |
---|
| 998 | pplay,pplev,zzlay,zzlev,z0, & |
---|
| 999 | pu,pv,zh,pq,tsurf,emis,qsurf, & |
---|
| 1000 | zdh,pdq,zflubid, & |
---|
| 1001 | zdudif,zdvdif,zdhdif,zdtsdif, & |
---|
[1524] | 1002 | sensibFlux,q2,zdqdif,zdqsdif, & |
---|
[1477] | 1003 | taux,tauy,lastcall) |
---|
[253] | 1004 | |
---|
[1477] | 1005 | zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only |
---|
[1524] | 1006 | zdqevap(1:ngrid,1:nlayer)=0. |
---|
[594] | 1007 | |
---|
[1477] | 1008 | end if !end of 'UseTurbDiff' |
---|
[594] | 1009 | |
---|
[1836] | 1010 | !!! this is always done, except for turbulence-resolving simulations |
---|
| 1011 | if (.not. turb_resolved) then |
---|
| 1012 | pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer) |
---|
| 1013 | pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer) |
---|
| 1014 | pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer) |
---|
| 1015 | zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid) |
---|
| 1016 | endif |
---|
[1477] | 1017 | |
---|
[1297] | 1018 | if(ok_slab_ocean)then |
---|
| 1019 | flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid)) |
---|
| 1020 | endif |
---|
| 1021 | |
---|
| 1022 | |
---|
[253] | 1023 | if (tracer) then |
---|
[1308] | 1024 | pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq) |
---|
[787] | 1025 | dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq) |
---|
[253] | 1026 | end if ! of if (tracer) |
---|
| 1027 | |
---|
[1477] | 1028 | |
---|
| 1029 | ! test energy conservation |
---|
[253] | 1030 | !------------------------- |
---|
| 1031 | if(enertest)then |
---|
[1477] | 1032 | |
---|
[1524] | 1033 | dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:) |
---|
[253] | 1034 | do ig = 1, ngrid |
---|
[1524] | 1035 | dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground |
---|
| 1036 | dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground |
---|
[253] | 1037 | enddo |
---|
[1477] | 1038 | |
---|
[1542] | 1039 | call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot) |
---|
[1524] | 1040 | dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:) |
---|
[1542] | 1041 | call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots) |
---|
| 1042 | call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux) |
---|
[1477] | 1043 | |
---|
[1524] | 1044 | if (is_master) then |
---|
[1477] | 1045 | |
---|
| 1046 | if (UseTurbDiff) then |
---|
[1524] | 1047 | print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2' |
---|
| 1048 | print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2' |
---|
[1477] | 1049 | print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2' |
---|
[1524] | 1050 | else |
---|
| 1051 | print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2' |
---|
| 1052 | print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2' |
---|
| 1053 | print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2' |
---|
| 1054 | end if |
---|
| 1055 | endif ! end of 'is_master' |
---|
[1477] | 1056 | |
---|
| 1057 | ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere. |
---|
| 1058 | endif ! end of 'enertest' |
---|
[253] | 1059 | |
---|
[1477] | 1060 | |
---|
| 1061 | ! Test water conservation. |
---|
[253] | 1062 | if(watertest.and.water)then |
---|
[1477] | 1063 | |
---|
[1524] | 1064 | call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp) |
---|
[1542] | 1065 | call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*cell_area(:)*ptimestep/totarea_planet,dWtots_tmp) |
---|
[253] | 1066 | do ig = 1, ngrid |
---|
[1524] | 1067 | vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap)) |
---|
| 1068 | enddo |
---|
| 1069 | call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) |
---|
[1542] | 1070 | call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)*cell_area(:)*ptimestep/totarea_planet,dWtots) |
---|
[1524] | 1071 | dWtot = dWtot + dWtot_tmp |
---|
| 1072 | dWtots = dWtots + dWtots_tmp |
---|
[651] | 1073 | do ig = 1, ngrid |
---|
[1524] | 1074 | vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice)) |
---|
| 1075 | enddo |
---|
| 1076 | call planetwide_maxval(vdifcncons(:),nconsMAX) |
---|
[253] | 1077 | |
---|
[1524] | 1078 | if (is_master) then |
---|
[1477] | 1079 | print*,'---------------------------------------------------------------' |
---|
| 1080 | print*,'In difv atmospheric water change =',dWtot,' kg m-2' |
---|
| 1081 | print*,'In difv surface water change =',dWtots,' kg m-2' |
---|
| 1082 | print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2' |
---|
| 1083 | print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1' |
---|
[1524] | 1084 | endif |
---|
[253] | 1085 | |
---|
[1477] | 1086 | endif ! end of 'watertest' |
---|
[253] | 1087 | !------------------------- |
---|
| 1088 | |
---|
[1477] | 1089 | else ! calldifv |
---|
[253] | 1090 | |
---|
| 1091 | if(.not.newtonian)then |
---|
| 1092 | |
---|
[787] | 1093 | zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid) |
---|
[253] | 1094 | |
---|
| 1095 | endif |
---|
| 1096 | |
---|
[1477] | 1097 | endif ! end of 'calldifv' |
---|
[253] | 1098 | |
---|
| 1099 | |
---|
[1477] | 1100 | !---------------------------------- |
---|
| 1101 | ! IV. Dry convective adjustment : |
---|
| 1102 | !---------------------------------- |
---|
[253] | 1103 | |
---|
| 1104 | if(calladj) then |
---|
| 1105 | |
---|
[1308] | 1106 | zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer) |
---|
| 1107 | zduadj(1:ngrid,1:nlayer)=0.0 |
---|
| 1108 | zdvadj(1:ngrid,1:nlayer)=0.0 |
---|
| 1109 | zdhadj(1:ngrid,1:nlayer)=0.0 |
---|
[253] | 1110 | |
---|
| 1111 | |
---|
[1477] | 1112 | call convadj(ngrid,nlayer,nq,ptimestep, & |
---|
| 1113 | pplay,pplev,zpopsk, & |
---|
| 1114 | pu,pv,zh,pq, & |
---|
| 1115 | pdu,pdv,zdh,pdq, & |
---|
| 1116 | zduadj,zdvadj,zdhadj, & |
---|
| 1117 | zdqadj) |
---|
[253] | 1118 | |
---|
[1308] | 1119 | pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer) |
---|
| 1120 | pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer) |
---|
| 1121 | pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) |
---|
| 1122 | zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only |
---|
[1283] | 1123 | |
---|
[253] | 1124 | if(tracer) then |
---|
[1308] | 1125 | pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq) |
---|
[253] | 1126 | end if |
---|
| 1127 | |
---|
[1477] | 1128 | ! Test energy conservation |
---|
[253] | 1129 | if(enertest)then |
---|
[1524] | 1130 | call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot) |
---|
[1295] | 1131 | if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2' |
---|
[253] | 1132 | endif |
---|
| 1133 | |
---|
[1477] | 1134 | ! Test water conservation |
---|
[253] | 1135 | if(watertest)then |
---|
[1524] | 1136 | call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp) |
---|
[253] | 1137 | do ig = 1, ngrid |
---|
[1524] | 1138 | cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap)) |
---|
| 1139 | enddo |
---|
| 1140 | call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) |
---|
| 1141 | dWtot = dWtot + dWtot_tmp |
---|
[651] | 1142 | do ig = 1, ngrid |
---|
[1524] | 1143 | cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice)) |
---|
| 1144 | enddo |
---|
| 1145 | call planetwide_maxval(cadjncons(:),nconsMAX) |
---|
[253] | 1146 | |
---|
[1295] | 1147 | if (is_master) then |
---|
[1524] | 1148 | print*,'In convadj atmospheric water change =',dWtot,' kg m-2' |
---|
[1477] | 1149 | print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1' |
---|
[1524] | 1150 | endif |
---|
[1477] | 1151 | |
---|
| 1152 | endif ! end of 'watertest' |
---|
[787] | 1153 | |
---|
[1477] | 1154 | endif ! end of 'calladj' |
---|
[253] | 1155 | |
---|
[1477] | 1156 | !----------------------------------------------- |
---|
| 1157 | ! V. Carbon dioxide condensation-sublimation : |
---|
| 1158 | !----------------------------------------------- |
---|
[253] | 1159 | |
---|
| 1160 | if (co2cond) then |
---|
[1477] | 1161 | |
---|
[253] | 1162 | if (.not.tracer) then |
---|
| 1163 | print*,'We need a CO2 ice tracer to condense CO2' |
---|
| 1164 | call abort |
---|
| 1165 | endif |
---|
[1477] | 1166 | call condense_co2(ngrid,nlayer,nq,ptimestep, & |
---|
| 1167 | capcal,pplay,pplev,tsurf,pt, & |
---|
[1485] | 1168 | pdt,zdtsurf,pq,pdq, & |
---|
| 1169 | qsurf,zdqsurfc,albedo,emis, & |
---|
[1482] | 1170 | albedo_bareground,albedo_co2_ice_SPECTV, & |
---|
[1485] | 1171 | zdtc,zdtsurfc,pdpsrf,zdqc) |
---|
[253] | 1172 | |
---|
[1484] | 1173 | pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer) |
---|
| 1174 | zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid) |
---|
[728] | 1175 | |
---|
[1484] | 1176 | pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq) |
---|
| 1177 | dqsurf(1:ngrid,igcm_co2_ice) = dqsurf(1:ngrid,igcm_co2_ice) + zdqsurfc(1:ngrid) |
---|
[253] | 1178 | |
---|
| 1179 | ! test energy conservation |
---|
| 1180 | if(enertest)then |
---|
[1524] | 1181 | call planetwide_sumval(cpp*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot) |
---|
[1542] | 1182 | call planetwide_sumval(capcal(:)*zdtsurfc(:)*cell_area(:)/totarea_planet,dEtots) |
---|
[1524] | 1183 | if (is_master) then |
---|
| 1184 | print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2' |
---|
[1477] | 1185 | print*,'In co2cloud surface energy change =',dEtots,' W m-2' |
---|
[1524] | 1186 | endif |
---|
[253] | 1187 | endif |
---|
| 1188 | |
---|
[1477] | 1189 | endif ! end of 'co2cond' |
---|
[253] | 1190 | |
---|
| 1191 | |
---|
[1477] | 1192 | !--------------------------------------------- |
---|
| 1193 | ! VI. Specific parameterizations for tracers |
---|
| 1194 | !--------------------------------------------- |
---|
[253] | 1195 | |
---|
[1477] | 1196 | if (tracer) then |
---|
| 1197 | |
---|
| 1198 | ! --------------------- |
---|
| 1199 | ! VI.1. Water and ice |
---|
| 1200 | ! --------------------- |
---|
[253] | 1201 | if (water) then |
---|
| 1202 | |
---|
[1477] | 1203 | ! Water ice condensation in the atmosphere |
---|
[728] | 1204 | if(watercond.and.(RLVTT.gt.1.e-8))then |
---|
[253] | 1205 | |
---|
[1524] | 1206 | dqmoist(1:ngrid,1:nlayer,1:nq)=0. |
---|
| 1207 | dtmoist(1:ngrid,1:nlayer)=0. |
---|
[1477] | 1208 | |
---|
| 1209 | ! Moist Convective Adjustment. |
---|
| 1210 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1524] | 1211 | call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man) |
---|
[253] | 1212 | |
---|
[1477] | 1213 | pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) & |
---|
[1524] | 1214 | + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) |
---|
[1477] | 1215 | pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) & |
---|
| 1216 | + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
[1308] | 1217 | pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer) |
---|
[728] | 1218 | |
---|
[1477] | 1219 | ! Test energy conservation. |
---|
[253] | 1220 | if(enertest)then |
---|
[1477] | 1221 | |
---|
[1524] | 1222 | call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot) |
---|
| 1223 | call planetwide_maxval(dtmoist(:,:),dtmoist_max) |
---|
| 1224 | call planetwide_minval(dtmoist(:,:),dtmoist_min) |
---|
| 1225 | madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:) |
---|
[787] | 1226 | do ig=1,ngrid |
---|
[728] | 1227 | madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:)) |
---|
[253] | 1228 | enddo |
---|
[1477] | 1229 | |
---|
[1524] | 1230 | if (is_master) then |
---|
[1477] | 1231 | print*,'In moistadj atmospheric energy change =',dEtot,' W m-2' |
---|
| 1232 | print*,'In moistadj MAX atmospheric energy change =',dtmoist_max*ptimestep,'K/step' |
---|
| 1233 | print*,'In moistadj MIN atmospheric energy change =',dtmoist_min*ptimestep,'K/step' |
---|
[1524] | 1234 | endif |
---|
[1477] | 1235 | |
---|
[1524] | 1236 | call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+ & |
---|
| 1237 | massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) |
---|
| 1238 | if (is_master) print*,'In moistadj atmospheric water change =',dWtot,' kg m-2' |
---|
| 1239 | |
---|
[1477] | 1240 | endif ! end of 'enertest' |
---|
[1524] | 1241 | |
---|
[253] | 1242 | |
---|
[1477] | 1243 | ! Large scale condensation/evaporation. |
---|
| 1244 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1308] | 1245 | call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc) |
---|
[253] | 1246 | |
---|
[1308] | 1247 | pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtlscale(1:ngrid,1:nlayer) |
---|
| 1248 | pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayer) |
---|
| 1249 | pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer) |
---|
[253] | 1250 | |
---|
[1477] | 1251 | ! Test energy conservation. |
---|
[253] | 1252 | if(enertest)then |
---|
[1016] | 1253 | lscaledEz(:,:) = cpp*mass(:,:)*dtlscale(:,:) |
---|
[787] | 1254 | do ig=1,ngrid |
---|
[728] | 1255 | lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:)) |
---|
[253] | 1256 | enddo |
---|
[1524] | 1257 | call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot) |
---|
[1477] | 1258 | |
---|
[1524] | 1259 | if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2' |
---|
[728] | 1260 | |
---|
[1477] | 1261 | ! Test water conservation. |
---|
[1524] | 1262 | call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+ & |
---|
| 1263 | SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot) |
---|
[1477] | 1264 | |
---|
[1524] | 1265 | if (is_master) print*,'In largescale atmospheric water change =',dWtot,' kg m-2' |
---|
[1477] | 1266 | endif ! end of 'enertest' |
---|
[253] | 1267 | |
---|
[1477] | 1268 | ! Compute cloud fraction. |
---|
[253] | 1269 | do l = 1, nlayer |
---|
[787] | 1270 | do ig=1,ngrid |
---|
[253] | 1271 | cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l)) |
---|
| 1272 | enddo |
---|
| 1273 | enddo |
---|
| 1274 | |
---|
[1477] | 1275 | endif ! end of 'watercond' |
---|
[253] | 1276 | |
---|
[1477] | 1277 | ! Water ice / liquid precipitation. |
---|
| 1278 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[728] | 1279 | if(waterrain)then |
---|
[253] | 1280 | |
---|
[1308] | 1281 | zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0 |
---|
[787] | 1282 | zdqsrain(1:ngrid) = 0.0 |
---|
| 1283 | zdqssnow(1:ngrid) = 0.0 |
---|
[253] | 1284 | |
---|
[1309] | 1285 | call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq, & |
---|
[1477] | 1286 | zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac) |
---|
[253] | 1287 | |
---|
[1308] | 1288 | pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) & |
---|
[1524] | 1289 | + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) |
---|
[1308] | 1290 | pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) & |
---|
[1524] | 1291 | + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
[1308] | 1292 | pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer) |
---|
[1477] | 1293 | |
---|
| 1294 | dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid) |
---|
| 1295 | dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid) |
---|
[253] | 1296 | |
---|
[1477] | 1297 | ! Test energy conservation. |
---|
[651] | 1298 | if(enertest)then |
---|
[1477] | 1299 | |
---|
[1524] | 1300 | call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot) |
---|
| 1301 | if (is_master) print*,'In rain atmospheric T energy change =',dEtot,' W m-2' |
---|
| 1302 | call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)/totarea_planet*RLVTT/cpp,dItot_tmp) |
---|
[1542] | 1303 | call planetwide_sumval(cell_area(:)*zdqssnow(:)/totarea_planet*RLVTT/cpp,dItot) |
---|
[1524] | 1304 | dItot = dItot + dItot_tmp |
---|
| 1305 | call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dVtot_tmp) |
---|
[1542] | 1306 | call planetwide_sumval(cell_area(:)*zdqsrain(:)/totarea_planet*RLVTT/cpp,dVtot) |
---|
[1524] | 1307 | dVtot = dVtot + dVtot_tmp |
---|
| 1308 | dEtot = dItot + dVtot |
---|
[1477] | 1309 | |
---|
[1524] | 1310 | if (is_master) then |
---|
[1477] | 1311 | print*,'In rain dItot =',dItot,' W m-2' |
---|
| 1312 | print*,'In rain dVtot =',dVtot,' W m-2' |
---|
| 1313 | print*,'In rain atmospheric L energy change =',dEtot,' W m-2' |
---|
[1524] | 1314 | endif |
---|
[253] | 1315 | |
---|
[1477] | 1316 | ! Test water conservation |
---|
[1524] | 1317 | call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+ & |
---|
| 1318 | massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) |
---|
[1542] | 1319 | call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*cell_area(:)*ptimestep/totarea_planet,dWtots) |
---|
[1477] | 1320 | |
---|
[1524] | 1321 | if (is_master) then |
---|
| 1322 | print*,'In rain atmospheric water change =',dWtot,' kg m-2' |
---|
| 1323 | print*,'In rain surface water change =',dWtots,' kg m-2' |
---|
| 1324 | print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2' |
---|
| 1325 | endif |
---|
[1477] | 1326 | |
---|
| 1327 | endif ! end of 'enertest' |
---|
[253] | 1328 | |
---|
[1477] | 1329 | end if ! enf of 'waterrain' |
---|
| 1330 | |
---|
| 1331 | end if ! end of 'water' |
---|
[253] | 1332 | |
---|
[1801] | 1333 | ! ------------------------- |
---|
| 1334 | ! VI.2. Photochemistry |
---|
[1477] | 1335 | ! ------------------------- |
---|
[1801] | 1336 | |
---|
| 1337 | IF (photochem) then |
---|
| 1338 | |
---|
| 1339 | DO ig=1,ngrid |
---|
| 1340 | array_zero1(ig)=0.0 |
---|
| 1341 | DO l=1,nlayer |
---|
| 1342 | array_zero2(ig,l)=0. |
---|
| 1343 | ENDDO |
---|
| 1344 | ENDDO |
---|
| 1345 | |
---|
| 1346 | call calchim_asis(ngrid,nlayer,nq, & |
---|
| 1347 | ptimestep,pplay,pplev,pt,pdt,dist_star,mu0, & |
---|
| 1348 | fract,zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, & |
---|
| 1349 | array_zero1,array_zero1, & |
---|
| 1350 | pu,pdu,pv,pdv,array_zero2,array_zero2) |
---|
| 1351 | |
---|
| 1352 | ! increment values of tracers: |
---|
| 1353 | DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry |
---|
| 1354 | ! tracers is zero anyways |
---|
| 1355 | DO l=1,nlayer |
---|
| 1356 | DO ig=1,ngrid |
---|
| 1357 | pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) |
---|
| 1358 | ENDDO |
---|
| 1359 | ENDDO |
---|
| 1360 | ENDDO ! of DO iq=1,nq |
---|
| 1361 | |
---|
| 1362 | |
---|
| 1363 | ! increment surface values of tracers: |
---|
| 1364 | DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry |
---|
| 1365 | ! tracers is zero anyways |
---|
| 1366 | DO ig=1,ngrid |
---|
| 1367 | ! dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) |
---|
| 1368 | ENDDO |
---|
| 1369 | ENDDO ! of DO iq=1,nq |
---|
| 1370 | |
---|
| 1371 | END IF ! of IF (photochem) |
---|
| 1372 | |
---|
| 1373 | |
---|
| 1374 | |
---|
[1477] | 1375 | ! ------------------------- |
---|
[1801] | 1376 | ! VI.3. Aerosol particles |
---|
| 1377 | ! ------------------------- |
---|
[253] | 1378 | |
---|
[1477] | 1379 | ! Sedimentation. |
---|
| 1380 | if (sedimentation) then |
---|
| 1381 | |
---|
| 1382 | zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0 |
---|
| 1383 | zdqssed(1:ngrid,1:nq) = 0.0 |
---|
[253] | 1384 | |
---|
[1477] | 1385 | if(watertest)then |
---|
| 1386 | |
---|
| 1387 | iq=igcm_h2o_ice |
---|
[1524] | 1388 | call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot) |
---|
| 1389 | call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots) |
---|
| 1390 | if (is_master) then |
---|
| 1391 | print*,'Before sedim pq =',dWtot,' kg m-2' |
---|
[1477] | 1392 | print*,'Before sedim pdq =',dWtots,' kg m-2' |
---|
[1524] | 1393 | endif |
---|
[1477] | 1394 | endif |
---|
| 1395 | |
---|
| 1396 | call callsedim(ngrid,nlayer,ptimestep, & |
---|
| 1397 | pplev,zzlev,pt,pdt,pq,pdq, & |
---|
| 1398 | zdqsed,zdqssed,nq) |
---|
[253] | 1399 | |
---|
[1477] | 1400 | if(watertest)then |
---|
| 1401 | iq=igcm_h2o_ice |
---|
[1524] | 1402 | call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot) |
---|
| 1403 | call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots) |
---|
| 1404 | if (is_master) then |
---|
| 1405 | print*,'After sedim pq =',dWtot,' kg m-2' |
---|
| 1406 | print*,'After sedim pdq =',dWtots,' kg m-2' |
---|
| 1407 | endif |
---|
[1477] | 1408 | endif |
---|
[253] | 1409 | |
---|
[1477] | 1410 | ! Whether it falls as rain or snow depends only on the surface temperature |
---|
| 1411 | pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq) |
---|
| 1412 | dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq) |
---|
[253] | 1413 | |
---|
[1477] | 1414 | ! Test water conservation |
---|
| 1415 | if(watertest)then |
---|
[1524] | 1416 | call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot) |
---|
[1542] | 1417 | call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*cell_area(:)*ptimestep/totarea_planet,dWtots) |
---|
[1524] | 1418 | if (is_master) then |
---|
| 1419 | print*,'In sedim atmospheric ice change =',dWtot,' kg m-2' |
---|
| 1420 | print*,'In sedim surface ice change =',dWtots,' kg m-2' |
---|
| 1421 | print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2' |
---|
| 1422 | endif |
---|
[1477] | 1423 | endif |
---|
[253] | 1424 | |
---|
[1477] | 1425 | end if ! end of 'sedimentation' |
---|
[253] | 1426 | |
---|
| 1427 | |
---|
[1477] | 1428 | ! --------------- |
---|
[1801] | 1429 | ! VI.4. Updates |
---|
[1477] | 1430 | ! --------------- |
---|
[253] | 1431 | |
---|
[1477] | 1432 | ! Updating Atmospheric Mass and Tracers budgets. |
---|
[728] | 1433 | if(mass_redistrib) then |
---|
| 1434 | |
---|
[1477] | 1435 | zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * & |
---|
[1524] | 1436 | ( zdqevap(1:ngrid,1:nlayer) & |
---|
| 1437 | + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) & |
---|
| 1438 | + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) & |
---|
| 1439 | + dqvaplscale(1:ngrid,1:nlayer) ) |
---|
[863] | 1440 | |
---|
| 1441 | do ig = 1, ngrid |
---|
[1524] | 1442 | zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer)) |
---|
[863] | 1443 | enddo |
---|
[728] | 1444 | |
---|
[1524] | 1445 | call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr) |
---|
| 1446 | call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col) |
---|
| 1447 | call writediagfi(ngrid,"mass","mass","kg/m2",3,mass) |
---|
[728] | 1448 | |
---|
[1524] | 1449 | call mass_redistribution(ngrid,nlayer,nq,ptimestep, & |
---|
[1477] | 1450 | rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf, & |
---|
[1524] | 1451 | pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, & |
---|
| 1452 | zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr) |
---|
| 1453 | |
---|
[1308] | 1454 | pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq) |
---|
[1477] | 1455 | dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq) |
---|
| 1456 | pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer) |
---|
| 1457 | pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer) |
---|
| 1458 | pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer) |
---|
[1524] | 1459 | pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid) |
---|
[1477] | 1460 | zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid) |
---|
[1524] | 1461 | |
---|
| 1462 | endif |
---|
[728] | 1463 | |
---|
[1477] | 1464 | ! ------------------ |
---|
[1801] | 1465 | ! VI.5. Slab Ocean |
---|
[1477] | 1466 | ! ------------------ |
---|
[728] | 1467 | |
---|
[1477] | 1468 | if (ok_slab_ocean)then |
---|
[728] | 1469 | |
---|
[1477] | 1470 | do ig=1,ngrid |
---|
| 1471 | qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice) |
---|
| 1472 | enddo |
---|
[1297] | 1473 | |
---|
[1477] | 1474 | call ocean_slab_ice(ptimestep, & |
---|
| 1475 | ngrid, knindex, tsea_ice, fluxrad, & |
---|
| 1476 | zdqssnow, qsurf(:,igcm_h2o_ice), & |
---|
| 1477 | - zdqsdif(:,igcm_h2o_vap), & |
---|
| 1478 | flux_sens_lat,tsea_ice, pctsrf_sic, & |
---|
| 1479 | taux,tauy,icount) |
---|
[1297] | 1480 | |
---|
| 1481 | |
---|
[1477] | 1482 | call ocean_slab_get_vars(ngrid,tslab, & |
---|
| 1483 | sea_ice, flux_o, & |
---|
| 1484 | flux_g, dt_hdiff, & |
---|
| 1485 | dt_ekman) |
---|
| 1486 | |
---|
[1297] | 1487 | do ig=1,ngrid |
---|
| 1488 | if (nint(rnat(ig)).eq.1)then |
---|
[1477] | 1489 | tslab(ig,1) = 0. |
---|
| 1490 | tslab(ig,2) = 0. |
---|
| 1491 | tsea_ice(ig) = 0. |
---|
| 1492 | sea_ice(ig) = 0. |
---|
| 1493 | pctsrf_sic(ig) = 0. |
---|
| 1494 | qsurf(ig,igcm_h2o_ice) = qsurfint(ig,igcm_h2o_ice) |
---|
[1297] | 1495 | endif |
---|
| 1496 | enddo |
---|
| 1497 | |
---|
[1477] | 1498 | endif ! end of 'ok_slab_ocean' |
---|
[1297] | 1499 | |
---|
[1477] | 1500 | ! ----------------------------- |
---|
[1801] | 1501 | ! VI.6. Surface Tracer Update |
---|
[1477] | 1502 | ! ----------------------------- |
---|
[1297] | 1503 | |
---|
[253] | 1504 | if(hydrology)then |
---|
[1297] | 1505 | |
---|
[1482] | 1506 | call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, & |
---|
| 1507 | capcal,albedo,albedo_bareground, & |
---|
[1524] | 1508 | albedo_snow_SPECTV,albedo_co2_ice_SPECTV, & |
---|
[1482] | 1509 | mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, & |
---|
[1477] | 1510 | sea_ice) |
---|
[253] | 1511 | |
---|
[1484] | 1512 | zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid) |
---|
| 1513 | dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + dqs_hyd(1:ngrid,1:nq) |
---|
| 1514 | |
---|
| 1515 | qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq) |
---|
[253] | 1516 | |
---|
[1477] | 1517 | ! Test energy conservation |
---|
[253] | 1518 | if(enertest)then |
---|
[1542] | 1519 | call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf_hyd(:)/totarea_planet,dEtots) |
---|
[1524] | 1520 | if (is_master) print*,'In hydrol surface energy change =',dEtots,' W m-2' |
---|
[253] | 1521 | endif |
---|
| 1522 | |
---|
| 1523 | ! test water conservation |
---|
| 1524 | if(watertest)then |
---|
[1542] | 1525 | call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)*cell_area(:)*ptimestep/totarea_planet,dWtots) |
---|
[1524] | 1526 | if (is_master) print*,'In hydrol surface ice change =',dWtots,' kg m-2' |
---|
[1542] | 1527 | call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*cell_area(:)*ptimestep/totarea_planet,dWtots) |
---|
[1524] | 1528 | if (is_master) then |
---|
[1477] | 1529 | print*,'In hydrol surface water change =',dWtots,' kg m-2' |
---|
| 1530 | print*,'---------------------------------------------------------------' |
---|
[1524] | 1531 | endif |
---|
[253] | 1532 | endif |
---|
| 1533 | |
---|
[1477] | 1534 | else ! of if (hydrology) |
---|
[253] | 1535 | |
---|
[1484] | 1536 | qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq) |
---|
[253] | 1537 | |
---|
[1477] | 1538 | end if ! of if (hydrology) |
---|
[253] | 1539 | |
---|
[1477] | 1540 | ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water |
---|
| 1541 | ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain. |
---|
[622] | 1542 | qsurf_hist(:,:) = qsurf(:,:) |
---|
[253] | 1543 | |
---|
| 1544 | if(ice_update)then |
---|
[787] | 1545 | ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice)) |
---|
[253] | 1546 | endif |
---|
| 1547 | |
---|
[1477] | 1548 | endif! end of if 'tracer' |
---|
[253] | 1549 | |
---|
| 1550 | |
---|
[1477] | 1551 | !------------------------------------------------ |
---|
| 1552 | ! VII. Surface and sub-surface soil temperature |
---|
| 1553 | !------------------------------------------------ |
---|
[253] | 1554 | |
---|
[1477] | 1555 | |
---|
| 1556 | ! Increment surface temperature |
---|
[1297] | 1557 | if(ok_slab_ocean)then |
---|
| 1558 | do ig=1,ngrid |
---|
| 1559 | if (nint(rnat(ig)).eq.0)then |
---|
| 1560 | zdtsurf(ig)= (tslab(ig,1) & |
---|
| 1561 | + pctsrf_sic(ig)*(tsea_ice(ig)-tslab(ig,1))-tsurf(ig))/ptimestep |
---|
| 1562 | endif |
---|
| 1563 | tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) |
---|
| 1564 | enddo |
---|
| 1565 | |
---|
| 1566 | else |
---|
| 1567 | tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) |
---|
[1477] | 1568 | endif ! end of 'ok_slab_ocean' |
---|
[1297] | 1569 | |
---|
[1477] | 1570 | |
---|
| 1571 | ! Compute soil temperatures and subsurface heat flux. |
---|
[253] | 1572 | if (callsoil) then |
---|
[787] | 1573 | call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, & |
---|
[1477] | 1574 | ptimestep,tsurf,tsoil,capcal,fluxgrd) |
---|
[253] | 1575 | endif |
---|
| 1576 | |
---|
[1297] | 1577 | |
---|
| 1578 | if (ok_slab_ocean) then |
---|
[1477] | 1579 | |
---|
| 1580 | do ig=1,ngrid |
---|
| 1581 | |
---|
| 1582 | fluxgrdocean(ig)=fluxgrd(ig) |
---|
| 1583 | if (nint(rnat(ig)).eq.0) then |
---|
[1297] | 1584 | capcal(ig)=capcalocean |
---|
| 1585 | fluxgrd(ig)=0. |
---|
| 1586 | fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1)) |
---|
[1477] | 1587 | do iq=1,nsoilmx |
---|
| 1588 | tsoil(ig,iq)=tsurf(ig) |
---|
| 1589 | enddo |
---|
| 1590 | if (pctsrf_sic(ig).gt.0.5) then |
---|
| 1591 | capcal(ig)=capcalseaice |
---|
| 1592 | if (qsurf(ig,igcm_h2o_ice).gt.0.) then |
---|
| 1593 | capcal(ig)=capcalsno |
---|
| 1594 | endif |
---|
| 1595 | endif |
---|
| 1596 | endif |
---|
| 1597 | |
---|
| 1598 | enddo |
---|
| 1599 | |
---|
| 1600 | endif !end of 'ok_slab_ocean' |
---|
[1297] | 1601 | |
---|
[1477] | 1602 | |
---|
| 1603 | ! Test energy conservation |
---|
[253] | 1604 | if(enertest)then |
---|
[1542] | 1605 | call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots) |
---|
[1524] | 1606 | if (is_master) print*,'Surface energy change =',dEtots,' W m-2' |
---|
[253] | 1607 | endif |
---|
| 1608 | |
---|
| 1609 | |
---|
[1477] | 1610 | !--------------------------------------------------- |
---|
| 1611 | ! VIII. Perform diagnostics and write output files |
---|
| 1612 | !--------------------------------------------------- |
---|
| 1613 | |
---|
| 1614 | ! Note : For output only: the actual model integration is performed in the dynamics. |
---|
| 1615 | |
---|
| 1616 | |
---|
[253] | 1617 | |
---|
[1477] | 1618 | ! Temperature, zonal and meridional winds. |
---|
[1308] | 1619 | zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep |
---|
| 1620 | zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep |
---|
| 1621 | zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep |
---|
[253] | 1622 | |
---|
[1477] | 1623 | ! Diagnostic. |
---|
[1637] | 1624 | zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep |
---|
[1308] | 1625 | ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer) |
---|
[253] | 1626 | |
---|
[1637] | 1627 | zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep |
---|
| 1628 | zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer) |
---|
| 1629 | |
---|
[253] | 1630 | if(firstcall)then |
---|
[1308] | 1631 | zdtdyn(1:ngrid,1:nlayer)=0.0 |
---|
[1637] | 1632 | zdudyn(1:ngrid,1:nlayer)=0.0 |
---|
[253] | 1633 | endif |
---|
| 1634 | |
---|
[1477] | 1635 | ! Dynamical heating diagnostic. |
---|
[253] | 1636 | do ig=1,ngrid |
---|
[1637] | 1637 | fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp |
---|
[253] | 1638 | enddo |
---|
| 1639 | |
---|
[1477] | 1640 | ! Tracers. |
---|
[1308] | 1641 | zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep |
---|
[253] | 1642 | |
---|
[1477] | 1643 | ! Surface pressure. |
---|
[787] | 1644 | ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep |
---|
[253] | 1645 | |
---|
| 1646 | |
---|
| 1647 | |
---|
[1477] | 1648 | ! Surface and soil temperature information |
---|
[1542] | 1649 | call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1) |
---|
[1295] | 1650 | call planetwide_minval(tsurf(:),Ts2) |
---|
| 1651 | call planetwide_maxval(tsurf(:),Ts3) |
---|
[253] | 1652 | if(callsoil)then |
---|
[1542] | 1653 | TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer |
---|
[1699] | 1654 | if (is_master) then |
---|
| 1655 | print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]' |
---|
| 1656 | print*,Ts1,Ts2,Ts3,TsS |
---|
| 1657 | end if |
---|
[959] | 1658 | else |
---|
[1699] | 1659 | if (is_master) then |
---|
| 1660 | print*,' ave[Tsurf] min[Tsurf] max[Tsurf]' |
---|
[1477] | 1661 | print*,Ts1,Ts2,Ts3 |
---|
[1524] | 1662 | endif |
---|
[959] | 1663 | end if |
---|
[253] | 1664 | |
---|
| 1665 | |
---|
[1477] | 1666 | ! Check the energy balance of the simulation during the run |
---|
[253] | 1667 | if(corrk)then |
---|
| 1668 | |
---|
[1542] | 1669 | call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR) |
---|
| 1670 | call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR) |
---|
| 1671 | call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR) |
---|
| 1672 | call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND) |
---|
| 1673 | call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN) |
---|
[787] | 1674 | do ig=1,ngrid |
---|
[253] | 1675 | if(fluxtop_dn(ig).lt.0.0)then |
---|
| 1676 | print*,'fluxtop_dn has gone crazy' |
---|
| 1677 | print*,'fluxtop_dn=',fluxtop_dn(ig) |
---|
| 1678 | print*,'tau_col=',tau_col(ig) |
---|
| 1679 | print*,'aerosol=',aerosol(ig,:,:) |
---|
| 1680 | print*,'temp= ',pt(ig,:) |
---|
| 1681 | print*,'pplay= ',pplay(ig,:) |
---|
| 1682 | call abort |
---|
| 1683 | endif |
---|
| 1684 | end do |
---|
| 1685 | |
---|
[787] | 1686 | if(ngrid.eq.1)then |
---|
[253] | 1687 | DYN=0.0 |
---|
| 1688 | endif |
---|
[1524] | 1689 | |
---|
| 1690 | if (is_master) then |
---|
[1477] | 1691 | print*,' ISR ASR OLR GND DYN [W m^-2]' |
---|
| 1692 | print*, ISR,ASR,OLR,GND,DYN |
---|
[1524] | 1693 | endif |
---|
[253] | 1694 | |
---|
[1295] | 1695 | if(enertest .and. is_master)then |
---|
[651] | 1696 | print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2' |
---|
| 1697 | print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2' |
---|
| 1698 | print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2' |
---|
[253] | 1699 | endif |
---|
| 1700 | |
---|
[1295] | 1701 | if(meanOLR .and. is_master)then |
---|
[1216] | 1702 | if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then |
---|
[253] | 1703 | ! to record global radiative balance |
---|
[588] | 1704 | open(92,file="rad_bal.out",form='formatted',position='append') |
---|
[651] | 1705 | write(92,*) zday,ISR,ASR,OLR |
---|
[253] | 1706 | close(92) |
---|
[588] | 1707 | open(93,file="tem_bal.out",form='formatted',position='append') |
---|
[1295] | 1708 | if(callsoil)then |
---|
[1524] | 1709 | write(93,*) zday,Ts1,Ts2,Ts3,TsS |
---|
| 1710 | else |
---|
| 1711 | write(93,*) zday,Ts1,Ts2,Ts3 |
---|
| 1712 | endif |
---|
[253] | 1713 | close(93) |
---|
| 1714 | endif |
---|
| 1715 | endif |
---|
| 1716 | |
---|
[1477] | 1717 | endif ! end of 'corrk' |
---|
[253] | 1718 | |
---|
[651] | 1719 | |
---|
[1477] | 1720 | ! Diagnostic to test radiative-convective timescales in code. |
---|
[253] | 1721 | if(testradtimes)then |
---|
[588] | 1722 | open(38,file="tau_phys.out",form='formatted',position='append') |
---|
[253] | 1723 | ig=1 |
---|
| 1724 | do l=1,nlayer |
---|
| 1725 | write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l) |
---|
| 1726 | enddo |
---|
| 1727 | close(38) |
---|
[726] | 1728 | print*,'As testradtimes enabled,' |
---|
| 1729 | print*,'exiting physics on first call' |
---|
[253] | 1730 | call abort |
---|
| 1731 | endif |
---|
| 1732 | |
---|
[1477] | 1733 | |
---|
| 1734 | ! Compute column amounts (kg m-2) if tracers are enabled. |
---|
[253] | 1735 | if(tracer)then |
---|
[787] | 1736 | qcol(1:ngrid,1:nq)=0.0 |
---|
[253] | 1737 | do iq=1,nq |
---|
[1477] | 1738 | do ig=1,ngrid |
---|
| 1739 | qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer)) |
---|
| 1740 | enddo |
---|
[253] | 1741 | enddo |
---|
| 1742 | |
---|
[1477] | 1743 | ! Generalised for arbitrary aerosols now. By LK |
---|
[787] | 1744 | reffcol(1:ngrid,1:naerkind)=0.0 |
---|
[728] | 1745 | if(co2cond.and.(iaero_co2.ne.0))then |
---|
[1308] | 1746 | call co2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_co2)) |
---|
[787] | 1747 | do ig=1,ngrid |
---|
[1308] | 1748 | reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayer,igcm_co2_ice)*reffrad(ig,1:nlayer,iaero_co2)*mass(ig,1:nlayer)) |
---|
[253] | 1749 | enddo |
---|
[728] | 1750 | endif |
---|
| 1751 | if(water.and.(iaero_h2o.ne.0))then |
---|
[1308] | 1752 | call h2o_reffrad(ngrid,nlayer,zq(1,1,igcm_h2o_ice),zt, & |
---|
[858] | 1753 | reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o)) |
---|
[787] | 1754 | do ig=1,ngrid |
---|
[1308] | 1755 | reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayer,igcm_h2o_ice)*reffrad(ig,1:nlayer,iaero_h2o)*mass(ig,1:nlayer)) |
---|
[728] | 1756 | enddo |
---|
| 1757 | endif |
---|
[253] | 1758 | |
---|
[1477] | 1759 | endif ! end of 'tracer' |
---|
[253] | 1760 | |
---|
| 1761 | |
---|
[1477] | 1762 | ! Test for water conservation. |
---|
[253] | 1763 | if(water)then |
---|
| 1764 | |
---|
[1542] | 1765 | call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf) |
---|
| 1766 | call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf) |
---|
| 1767 | call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol) |
---|
| 1768 | call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol) |
---|
[253] | 1769 | |
---|
[651] | 1770 | h2otot = icesrf + liqsrf + icecol + vapcol |
---|
[1524] | 1771 | |
---|
| 1772 | if (is_master) then |
---|
[1477] | 1773 | print*,' Total water amount [kg m^-2]: ',h2otot |
---|
| 1774 | print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] ' |
---|
| 1775 | print*, icesrf,liqsrf,icecol,vapcol |
---|
[1524] | 1776 | endif |
---|
[253] | 1777 | |
---|
[1295] | 1778 | if(meanOLR .and. is_master)then |
---|
[1216] | 1779 | if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then |
---|
[253] | 1780 | ! to record global water balance |
---|
[588] | 1781 | open(98,file="h2o_bal.out",form='formatted',position='append') |
---|
[651] | 1782 | write(98,*) zday,icesrf,liqsrf,icecol,vapcol |
---|
[253] | 1783 | close(98) |
---|
| 1784 | endif |
---|
| 1785 | endif |
---|
| 1786 | |
---|
| 1787 | endif |
---|
| 1788 | |
---|
| 1789 | |
---|
[1477] | 1790 | ! Calculate RH (Relative Humidity) for diagnostic. |
---|
[253] | 1791 | if(water)then |
---|
[1477] | 1792 | |
---|
[253] | 1793 | do l = 1, nlayer |
---|
[787] | 1794 | do ig=1,ngrid |
---|
[728] | 1795 | call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l)) |
---|
[253] | 1796 | RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l) |
---|
| 1797 | enddo |
---|
| 1798 | enddo |
---|
| 1799 | |
---|
[1477] | 1800 | ! Compute maximum possible H2O column amount (100% saturation). |
---|
[253] | 1801 | do ig=1,ngrid |
---|
[1477] | 1802 | H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:)) |
---|
[253] | 1803 | enddo |
---|
| 1804 | |
---|
[1477] | 1805 | endif ! end of 'water' |
---|
[253] | 1806 | |
---|
[996] | 1807 | |
---|
[1699] | 1808 | if (is_master) print*,'--> Ls =',zls*180./pi |
---|
[1477] | 1809 | |
---|
| 1810 | |
---|
| 1811 | !---------------------------------------------------------------------- |
---|
[253] | 1812 | ! Writing NetCDF file "RESTARTFI" at the end of the run |
---|
[1477] | 1813 | !---------------------------------------------------------------------- |
---|
| 1814 | |
---|
[253] | 1815 | ! Note: 'restartfi' is stored just before dynamics are stored |
---|
| 1816 | ! in 'restart'. Between now and the writting of 'restart', |
---|
| 1817 | ! there will have been the itau=itau+1 instruction and |
---|
| 1818 | ! a reset of 'time' (lastacll = .true. when itau+1= itaufin) |
---|
| 1819 | ! thus we store for time=time+dtvr |
---|
| 1820 | |
---|
| 1821 | |
---|
| 1822 | |
---|
[1477] | 1823 | if(lastcall) then |
---|
| 1824 | ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) |
---|
[305] | 1825 | |
---|
[1477] | 1826 | ! Update surface ice distribution to iterate to steady state if requested |
---|
| 1827 | if(ice_update)then |
---|
[253] | 1828 | |
---|
[1477] | 1829 | do ig=1,ngrid |
---|
[305] | 1830 | |
---|
[1477] | 1831 | delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig)) |
---|
| 1832 | |
---|
| 1833 | ! add multiple years of evolution |
---|
| 1834 | qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep |
---|
[305] | 1835 | |
---|
[1477] | 1836 | ! if ice has gone -ve, set to zero |
---|
| 1837 | if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then |
---|
| 1838 | qsurf_hist(ig,igcm_h2o_ice) = 0.0 |
---|
| 1839 | endif |
---|
[305] | 1840 | |
---|
[1477] | 1841 | ! if ice is seasonal, set to zero (NEW) |
---|
| 1842 | if(ice_min(ig).lt.0.01)then |
---|
| 1843 | qsurf_hist(ig,igcm_h2o_ice) = 0.0 |
---|
| 1844 | endif |
---|
[253] | 1845 | |
---|
[1477] | 1846 | enddo |
---|
[305] | 1847 | |
---|
[1477] | 1848 | ! enforce ice conservation |
---|
[1542] | 1849 | ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*cell_area(:) )/SUM(cell_area(:)) |
---|
[1477] | 1850 | qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot) |
---|
[305] | 1851 | |
---|
[1477] | 1852 | endif |
---|
[1836] | 1853 | #ifndef MESOSCALE |
---|
| 1854 | |
---|
[1477] | 1855 | if (ngrid.ne.1) then |
---|
| 1856 | write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin |
---|
| 1857 | |
---|
| 1858 | call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, & |
---|
| 1859 | ptimestep,ztime_fin, & |
---|
| 1860 | tsurf,tsoil,emis,q2,qsurf_hist, & |
---|
| 1861 | cloudfrac,totcloudfrac,hice, & |
---|
| 1862 | rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) |
---|
| 1863 | endif |
---|
[1836] | 1864 | #endif |
---|
[1477] | 1865 | if(ok_slab_ocean) then |
---|
| 1866 | call ocean_slab_final!(tslab, seaice) |
---|
| 1867 | end if |
---|
[1297] | 1868 | |
---|
[1682] | 1869 | endif ! end of 'lastcall' |
---|
[253] | 1870 | |
---|
[861] | 1871 | |
---|
[1477] | 1872 | !----------------------------------- |
---|
[253] | 1873 | ! Saving statistics : |
---|
[1477] | 1874 | !----------------------------------- |
---|
[253] | 1875 | |
---|
[1477] | 1876 | ! Note :("stats" stores and accumulates 8 key variables in file "stats.nc" |
---|
| 1877 | ! which can later be used to make the statistic files of the run: |
---|
| 1878 | ! "stats") only possible in 3D runs !!! |
---|
| 1879 | |
---|
[253] | 1880 | |
---|
[1477] | 1881 | if (callstats) then |
---|
[253] | 1882 | |
---|
[1477] | 1883 | call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) |
---|
| 1884 | call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf) |
---|
| 1885 | call wstats(ngrid,"fluxsurf_lw", & |
---|
| 1886 | "Thermal IR radiative flux to surface","W.m-2",2, & |
---|
| 1887 | fluxsurf_lw) |
---|
| 1888 | call wstats(ngrid,"fluxtop_lw", & |
---|
| 1889 | "Thermal IR radiative flux to space","W.m-2",2, & |
---|
| 1890 | fluxtop_lw) |
---|
| 1891 | |
---|
[253] | 1892 | ! call wstats(ngrid,"fluxsurf_sw", & |
---|
| 1893 | ! "Solar radiative flux to surface","W.m-2",2, & |
---|
[1477] | 1894 | ! fluxsurf_sw_tot) |
---|
[253] | 1895 | ! call wstats(ngrid,"fluxtop_sw", & |
---|
| 1896 | ! "Solar radiative flux to space","W.m-2",2, & |
---|
| 1897 | ! fluxtop_sw_tot) |
---|
[526] | 1898 | |
---|
[253] | 1899 | |
---|
[1477] | 1900 | call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn) |
---|
| 1901 | call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) |
---|
| 1902 | call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) |
---|
[1482] | 1903 | !call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) |
---|
| 1904 | !call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) |
---|
[1477] | 1905 | call wstats(ngrid,"p","Pressure","Pa",3,pplay) |
---|
| 1906 | call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) |
---|
| 1907 | call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) |
---|
| 1908 | call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv) |
---|
| 1909 | call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw) |
---|
| 1910 | call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2) |
---|
[526] | 1911 | |
---|
[1477] | 1912 | if (tracer) then |
---|
| 1913 | do iq=1,nq |
---|
| 1914 | call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) |
---|
| 1915 | call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & |
---|
| 1916 | 'kg m^-2',2,qsurf(1,iq) ) |
---|
| 1917 | call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & |
---|
[526] | 1918 | 'kg m^-2',2,qcol(1,iq) ) |
---|
[1477] | 1919 | |
---|
| 1920 | ! call wstats(ngrid,trim(noms(iq))//'_reff', & |
---|
[726] | 1921 | ! trim(noms(iq))//'_reff', & |
---|
| 1922 | ! 'm',3,reffrad(1,1,iq)) |
---|
[1477] | 1923 | |
---|
| 1924 | end do |
---|
| 1925 | |
---|
[253] | 1926 | if (water) then |
---|
[1308] | 1927 | vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap) |
---|
[1477] | 1928 | call wstats(ngrid,"vmr_h2ovapor", & |
---|
| 1929 | "H2O vapour volume mixing ratio","mol/mol", & |
---|
| 1930 | 3,vmr) |
---|
| 1931 | endif |
---|
[253] | 1932 | |
---|
[1477] | 1933 | endif ! end of 'tracer' |
---|
[253] | 1934 | |
---|
[1477] | 1935 | if(watercond.and.CLFvarying)then |
---|
| 1936 | call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man) |
---|
| 1937 | call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc) |
---|
| 1938 | call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac) |
---|
| 1939 | call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac) |
---|
| 1940 | call wstats(ngrid,"RH","relative humidity"," ",3,RH) |
---|
| 1941 | endif |
---|
[1297] | 1942 | |
---|
[1477] | 1943 | if (ok_slab_ocean) then |
---|
[1297] | 1944 | call wstats(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1)) |
---|
| 1945 | call wstats(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2)) |
---|
| 1946 | call wstats(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1)) |
---|
| 1947 | call wstats(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2)) |
---|
| 1948 | call wstats(ngrid,"tslab1","tslab1","K",2,tslab(:,1)) |
---|
| 1949 | call wstats(ngrid,"tslab2","tslab2","K",2,tslab(:,2)) |
---|
| 1950 | call wstats(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic) |
---|
| 1951 | call wstats(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice) |
---|
| 1952 | call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice) |
---|
| 1953 | call wstats(ngrid,"rnat","nature of the surface","",2,rnat) |
---|
[1477] | 1954 | endif |
---|
[1297] | 1955 | |
---|
[1477] | 1956 | if(lastcall) then |
---|
| 1957 | write (*,*) "Writing stats..." |
---|
| 1958 | call mkstats(ierr) |
---|
| 1959 | endif |
---|
| 1960 | |
---|
| 1961 | endif ! end of 'callstats' |
---|
[253] | 1962 | |
---|
[1836] | 1963 | #ifndef MESOSCALE |
---|
| 1964 | |
---|
[1477] | 1965 | !----------------------------------------------------------------------------------------------------- |
---|
| 1966 | ! OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic |
---|
| 1967 | ! |
---|
| 1968 | ! Note 1 : output with period "ecritphy", set in "run.def" |
---|
| 1969 | ! |
---|
| 1970 | ! Note 2 : writediagfi can also be called from any other subroutine for any variable, |
---|
| 1971 | ! but its preferable to keep all the calls in one place ... |
---|
| 1972 | !----------------------------------------------------------------------------------------------------- |
---|
[253] | 1973 | |
---|
[1477] | 1974 | call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi) |
---|
| 1975 | call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi) |
---|
| 1976 | call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi) |
---|
| 1977 | call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi) |
---|
| 1978 | call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf) |
---|
| 1979 | call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps) |
---|
| 1980 | call writediagfi(ngrid,"temp","temperature","K",3,zt) |
---|
| 1981 | call writediagfi(ngrid,"teta","potential temperature","K",3,zh) |
---|
| 1982 | call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu) |
---|
| 1983 | call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv) |
---|
| 1984 | call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw) |
---|
| 1985 | call writediagfi(ngrid,"p","Pressure","Pa",3,pplay) |
---|
| 1986 | |
---|
[965] | 1987 | ! Subsurface temperatures |
---|
[969] | 1988 | ! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf) |
---|
| 1989 | ! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil) |
---|
[965] | 1990 | |
---|
[1477] | 1991 | ! Oceanic layers |
---|
| 1992 | if(ok_slab_ocean) then |
---|
| 1993 | call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic) |
---|
| 1994 | call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice) |
---|
| 1995 | call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice) |
---|
| 1996 | call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1)) |
---|
| 1997 | call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2)) |
---|
| 1998 | call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1)) |
---|
| 1999 | call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2)) |
---|
| 2000 | call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1)) |
---|
| 2001 | call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2)) |
---|
| 2002 | call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat) |
---|
| 2003 | call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux) |
---|
| 2004 | call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT) |
---|
| 2005 | endif |
---|
| 2006 | |
---|
[1297] | 2007 | |
---|
[1477] | 2008 | ! Total energy balance diagnostics |
---|
| 2009 | if(callrad.and.(.not.newtonian))then |
---|
| 2010 | |
---|
[1482] | 2011 | !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) |
---|
| 2012 | !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) |
---|
[1477] | 2013 | call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn) |
---|
| 2014 | call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) |
---|
| 2015 | call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) |
---|
| 2016 | call writediagfi(ngrid,"shad","rings"," ", 2, fract) |
---|
| 2017 | |
---|
[1016] | 2018 | ! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1) |
---|
| 2019 | ! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1) |
---|
| 2020 | ! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw) |
---|
| 2021 | ! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw) |
---|
| 2022 | ! call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1) |
---|
| 2023 | ! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1) |
---|
[1477] | 2024 | |
---|
| 2025 | if(ok_slab_ocean) then |
---|
| 2026 | call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean) |
---|
| 2027 | else |
---|
| 2028 | call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd) |
---|
| 2029 | endif |
---|
| 2030 | |
---|
| 2031 | call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn) |
---|
| 2032 | |
---|
| 2033 | endif ! end of 'callrad' |
---|
[1524] | 2034 | |
---|
[1477] | 2035 | if(enertest) then |
---|
| 2036 | |
---|
[1524] | 2037 | if (calldifv) then |
---|
[1477] | 2038 | |
---|
[1524] | 2039 | call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2) |
---|
[1477] | 2040 | call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux) |
---|
| 2041 | |
---|
[1524] | 2042 | ! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff) |
---|
| 2043 | ! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff) |
---|
| 2044 | ! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs) |
---|
| 2045 | |
---|
| 2046 | endif |
---|
[1477] | 2047 | |
---|
[1524] | 2048 | if (corrk) then |
---|
| 2049 | call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw) |
---|
| 2050 | call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw) |
---|
| 2051 | endif |
---|
[1477] | 2052 | |
---|
| 2053 | if(watercond) then |
---|
| 2054 | |
---|
[1524] | 2055 | call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) |
---|
| 2056 | call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE) |
---|
| 2057 | call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat) |
---|
[1477] | 2058 | |
---|
[1524] | 2059 | ! call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) |
---|
| 2060 | ! call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz) |
---|
| 2061 | ! call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol) |
---|
[253] | 2062 | |
---|
[1477] | 2063 | endif |
---|
| 2064 | |
---|
| 2065 | endif ! end of 'enertest' |
---|
[253] | 2066 | |
---|
[1477] | 2067 | |
---|
| 2068 | ! Temporary inclusions for heating diagnostics. |
---|
| 2069 | !call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw) |
---|
| 2070 | !call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw) |
---|
| 2071 | !call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad) |
---|
| 2072 | ! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn) |
---|
| 2073 | |
---|
| 2074 | ! For Debugging. |
---|
[368] | 2075 | !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat)) |
---|
[253] | 2076 | !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi) |
---|
[1477] | 2077 | |
---|
[253] | 2078 | |
---|
[1477] | 2079 | ! Output aerosols. |
---|
| 2080 | if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) & |
---|
[1524] | 2081 | call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2)) |
---|
[1477] | 2082 | if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & |
---|
[1524] | 2083 | call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o)) |
---|
[1477] | 2084 | if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) & |
---|
[1524] | 2085 | call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2)) |
---|
[1477] | 2086 | if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & |
---|
[1524] | 2087 | call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o)) |
---|
[253] | 2088 | |
---|
[1477] | 2089 | ! Output tracers. |
---|
| 2090 | if (tracer) then |
---|
| 2091 | |
---|
| 2092 | do iq=1,nq |
---|
| 2093 | call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) |
---|
| 2094 | call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & |
---|
| 2095 | 'kg m^-2',2,qsurf_hist(1,iq) ) |
---|
| 2096 | call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & |
---|
| 2097 | 'kg m^-2',2,qcol(1,iq) ) |
---|
| 2098 | |
---|
[787] | 2099 | ! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & |
---|
[1477] | 2100 | ! 'kg m^-2',2,qsurf(1,iq) ) |
---|
[253] | 2101 | |
---|
[1477] | 2102 | if(watercond.or.CLFvarying)then |
---|
| 2103 | call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man) |
---|
| 2104 | call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc) |
---|
| 2105 | call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac) |
---|
| 2106 | call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac) |
---|
[1524] | 2107 | call writediagfi(ngrid,"RH","relative humidity"," ",3,RH) |
---|
[1477] | 2108 | endif |
---|
[253] | 2109 | |
---|
[1477] | 2110 | if(waterrain)then |
---|
| 2111 | call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain) |
---|
| 2112 | call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow) |
---|
| 2113 | endif |
---|
[253] | 2114 | |
---|
[1477] | 2115 | if((hydrology).and.(.not.ok_slab_ocean))then |
---|
| 2116 | call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice) |
---|
| 2117 | endif |
---|
[253] | 2118 | |
---|
[1477] | 2119 | if(ice_update)then |
---|
| 2120 | call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min) |
---|
| 2121 | call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial) |
---|
| 2122 | endif |
---|
[253] | 2123 | |
---|
[1477] | 2124 | do ig=1,ngrid |
---|
| 2125 | if(tau_col(ig).gt.1.e3)then |
---|
| 2126 | print*,'WARNING: tau_col=',tau_col(ig) |
---|
| 2127 | print*,'at ig=',ig,'in PHYSIQ' |
---|
| 2128 | endif |
---|
| 2129 | end do |
---|
[253] | 2130 | |
---|
[1477] | 2131 | call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col) |
---|
[253] | 2132 | |
---|
[1477] | 2133 | enddo ! end of 'nq' loop |
---|
| 2134 | |
---|
| 2135 | endif ! end of 'tracer' |
---|
[253] | 2136 | |
---|
[1477] | 2137 | |
---|
| 2138 | ! Output spectrum. |
---|
[526] | 2139 | if(specOLR.and.corrk)then |
---|
[728] | 2140 | call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu) |
---|
| 2141 | call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu) |
---|
[526] | 2142 | endif |
---|
[253] | 2143 | |
---|
[1836] | 2144 | #else |
---|
| 2145 | comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer) |
---|
| 2146 | comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer) |
---|
| 2147 | comm_CLOUDFRAC(1:ngrid,1:nlayer)=cloudfrac(1:ngrid,1:nlayer) |
---|
| 2148 | comm_TOTCLOUDFRAC(1:ngrid)=totcloudfrac(1:ngrid) |
---|
| 2149 | comm_RAIN(1:ngrid,1:nlayer)=zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) |
---|
| 2150 | comm_SNOW(1:ngrid,1:nlayer)=zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
| 2151 | comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid) |
---|
| 2152 | comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid) |
---|
| 2153 | comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid) |
---|
| 2154 | comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid) |
---|
| 2155 | comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid) |
---|
| 2156 | comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid) |
---|
| 2157 | comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid) |
---|
| 2158 | comm_LSCEZ(1:ngrid,1:nlayer)=lscaledEz(1:ngrid,1:nlayer) |
---|
| 2159 | comm_H2OICE_REFF(1:ngrid,1:nlayer)=reffrad(1:ngrid,1:nlayer,iaero_h2o) |
---|
| 2160 | sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ???? |
---|
| 2161 | #endif |
---|
| 2162 | |
---|
[1622] | 2163 | ! XIOS outputs |
---|
| 2164 | #ifdef CPP_XIOS |
---|
| 2165 | ! Send fields to XIOS: (NB these fields must also be defined as |
---|
| 2166 | ! <field id="..." /> in context_lmdz_physics.xml to be correctly used) |
---|
[1626] | 2167 | CALL send_xios_field("ls",zls) |
---|
| 2168 | |
---|
[1622] | 2169 | CALL send_xios_field("ps",ps) |
---|
| 2170 | CALL send_xios_field("area",cell_area) |
---|
| 2171 | |
---|
| 2172 | CALL send_xios_field("temperature",zt) |
---|
| 2173 | CALL send_xios_field("u",zu) |
---|
| 2174 | CALL send_xios_field("v",zv) |
---|
[1682] | 2175 | |
---|
| 2176 | if (lastcall.and.is_omp_master) then |
---|
| 2177 | write(*,*) "physiq: call xios_context_finalize" |
---|
| 2178 | call xios_context_finalize |
---|
| 2179 | endif |
---|
[1622] | 2180 | #endif |
---|
| 2181 | |
---|
[253] | 2182 | icount=icount+1 |
---|
| 2183 | |
---|
[751] | 2184 | end subroutine physiq |
---|
[1549] | 2185 | |
---|
| 2186 | end module physiq_mod |
---|