source: trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 @ 2757

Last change on this file since 2757 was 2736, checked in by aslmd, 3 years ago

Changed of callcorrk to physiq_mod

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