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