source: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 @ 3640

Last change on this file since 3640 was 3640, checked in by tbertrand, 4 months ago

Pluto: fix to take into account change in soil thermal inertia with time + fix for albedomap in newstart
TB

  • Property svn:executable set to *
File size: 98.8 KB
Line 
1      module physiq_mod
2
3      implicit none
4
5      contains
6
7      subroutine physiq(ngrid,nlayer,nq,   &
8                  firstcall,lastcall,      &
9                  pday,ptime,ptimestep,    &
10                  pplev,pplay,pphi,        &
11                  pu,pv,pt,pq,             &
12                  flxw,                    &
13                  pdu,pdv,pdt,pdq,pdpsrf)
14
15!!
16      use write_field_phy, only: Writefield_phy
17!!
18      use ioipsl_getin_p_mod, only: getin_p
19      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
20      use gases_h, only: gnom, gfrac
21      use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
22      use suaer_corrk_mod, only: suaer_corrk
23      use aerosol_mod, only: i_haze, haze_prof
24      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
25                           dryness
26      use comdiurn_h, only: coslat, sinlat, coslon, sinlon
27      use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
28      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
29      use geometry_mod, only: latitude, longitude, cell_area, &
30                          cell_area_for_lonlat_outputs
31      USE comgeomfi_h, only: totarea, totarea_planet
32      USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
33                          igcm_n2,igcm_ch4_gas,igcm_ch4_ice,igcm_haze,&
34                          igcm_co_gas,igcm_co_ice,igcm_prec_haze,lw_n2,lw_ch4,lw_co,&
35                          alpha_lift, alpha_devil, qextrhor, &
36                          nesp, is_chim
37      use time_phylmdz_mod, only: ecritphy, iphysiq, nday
38      use phyetat0_mod, only: phyetat0,tab_cntrl_mod
39      use wstats_mod, only: callstats, wstats, mkstats
40      use phyredem, only: physdem0, physdem1
41      use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
42      use mod_phys_lmdz_para, only : is_master
43      use planete_mod, only: apoastr, periastr, year_day, peri_day, &
44                            obliquit, z0, adjust, tpal
45      use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
46      use time_phylmdz_mod, only: daysec
47      use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv,        &
48                              callrad, callsoil, nosurf,                      &
49                              callconduct,callmolvis,callmoldiff,             &
50                              corrk, diagdtau,                                &
51                              diurnal, enertest, fat1au,                      &
52                              icetstep, intheat, iradia, kastprof,            &
53                              lwrite, mass_redistrib, meanOLR,                &
54                              fast,fasthaze,haze,metcloud,monoxcloud,         &
55                              n2cond,noseason_day,conservn2,conservch4,       &
56                              convergeps,kbo,triton,paleo,paleoyears,glaflow, &
57                              carbox, methane,condmetsurf,condcosurf,         &
58                              oldplutovdifc,oldplutocorrk,oldplutosedim,      &
59                              optichaze,haze_proffix,haze_radproffix,         &
60                              source_haze, tsurfmax, albmin_ch4,              &
61                              season, sedimentation,                          &
62                              specOLR,                                        &
63                              startphy_file, testradtimes,                    &
64                              tracer, UseTurbDiff,                            &
65                              global1d, szangle,                              &
66                              callmufi
67      use check_fields_mod, only: check_physics_fields
68      use conc_mod, only: rnew, cpnew, ini_conc_mod
69      use phys_state_var_mod
70      use callcorrk_mod, only: callcorrk
71      use callcorrk_pluto_mod, only: callcorrk_pluto
72      use vdifc_mod, only: vdifc
73      use vdifc_pluto_mod, only: vdifc_pluto
74      use turbdiff_mod, only: turbdiff
75      use turb_mod, only : q2,sensibFlux,turb_resolved
76      use mass_redistribution_mod, only: mass_redistribution
77      use datafile_mod, only: datadir
78      USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
79      use mod_phys_lmdz_omp_data, ONLY: is_omp_master
80      USE mod_grid_phy_lmdz, ONLY: regular_lonlat, grid_type, unstructured
81      ! Microphysical model (mp2m)
82      use mp2m_calmufi
83      use mp2m_diagnostics
84
85#ifdef CPP_XIOS
86      use xios_output_mod, only: initialize_xios_output, &
87                                 update_xios_timestep, &
88                                 send_xios_field
89      use wxios, only: wxios_context_init, xios_context_finalize
90#endif
91      USE mod_grid_phy_lmdz, ONLY: grid_type,unstructured,regular_lonlat
92      use write_output_mod, only: write_output
93
94      implicit none
95
96
97!==================================================================
98!
99!     Purpose
100!     -------
101!     Central subroutine for all the physics parameterisations in the
102!     Pluto model. Originally adapted from the Mars LMDZ model.
103!
104!     The model can be run with 1 (N2) or more tracer transport
105!     depending on the value of "tracer" in file "callphys.def".
106!
107!     It includes:
108!
109!      I. Initialization :
110!         I.1 Firstcall initializations.
111!         I.2 Initialization for every call to physiq.
112!
113!      II.1 Thermosphere
114!      II.2 Compute radiative transfer tendencies (longwave and shortwave) :
115!         II.2.a Option 1 : Call correlated-k radiative transfer scheme.
116!         II.2.b Option 2 : Atmosphere has no radiative effect.
117!
118!      III. Vertical diffusion (turbulent mixing)
119!
120!      IV. Convection :
121!         IV.a Dry convective adjusment
122!
123!      V. Condensation and sublimation of gases.
124!
125!      VI. Tracers
126!         VI.1. Microphysics / Aerosols and particles.
127!         VI.2. Updates (pressure variations, surface budget).
128!         VI.3. Surface Tracer Update.
129!
130!      VII. Surface and sub-surface soil temperature.
131!
132!      VIII. Perform diagnostics and write output files.
133!
134!
135!   arguments
136!   ---------
137!
138!   INPUT
139!   -----
140!
141!    ngrid                 Size of the horizontal grid.
142!    nlayer                Number of vertical layers.
143!    nq                    Number of advected fields.
144!
145!    firstcall             True at the first call.
146!    lastcall              True at the last call.
147!
148!    pday                  Number of days counted from the North. Spring equinoxe.
149!    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
150!    ptimestep             timestep (s).
151!
152!    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa).
153!    pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
154!    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2.s-2).
155!
156!    pu(ngrid,nlayer)      u, zonal component of the wind (ms-1).
157!    pv(ngrid,nlayer)      v, meridional component of the wind (ms-1).
158!
159!    pt(ngrid,nlayer)      Temperature (K).
160!
161!    pq(ngrid,nlayer,nq)   Advected fields.
162!
163!    pudyn(ngrid,nlayer)    \
164!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
165!    ptdyn(ngrid,nlayer)     / corresponding variables.
166!    pqdyn(ngrid,nlayer,nq) /
167!    flxw(ngrid,nlayer)      vertical mass flux (kg/s) at layer lower boundary
168!
169!   OUTPUT
170!   ------
171!
172!    pdu(ngrid,nlayer)        \
173!    pdv(ngrid,nlayer)         \  Temporal derivative of the corresponding
174!    pdt(ngrid,nlayer)         /  variables due to physical processes.
175!    pdq(ngrid,nlayer)        /
176!    pdpsrf(ngrid)             /
177!
178!
179!     Authors
180!     -------
181!           Frederic Hourdin        15/10/93
182!           Francois Forget        1994
183!           Christophe Hourdin        02/1997
184!           Subroutine completely rewritten by F. Forget (01/2000)
185!           Water ice clouds: Franck Montmessin (update 06/2003)
186!           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
187!           New correlated-k radiative scheme: R. Wordsworth (2009)
188!           Many specifically Martian subroutines removed: R. Wordsworth (2009)
189!           Improved water cycle: R. Wordsworth / B. Charnay (2010)
190!           To F90: R. Wordsworth (2010)
191!           New turbulent diffusion scheme: J. Leconte (2012)
192!           Loops converted to F90 matrix format: J. Leconte (2012)
193!           No more ngrid/nq, F90 commons and adaptation to parallel: A. Spiga (2012)
194!           Purge of the code : M. Turbet (2015)
195!           Photochemical core developped by F. Lefevre: B. Charnay (2017)
196!           Purge for Pluto model : A. Falco (2024)
197!           Adapting to Pluto : A. Falco, T. Bertrand (2024)
198!           Microphysical moment model: B. de Batz de Trenquelléon (2024)
199!==================================================================
200
201
202!    0.  Declarations :
203!    ------------------
204
205include "netcdf.inc"
206
207! Arguments :
208! -----------
209
210!   INPUTS:
211!   -------
212
213      integer,intent(in) :: ngrid             ! Number of atmospheric columns.
214      integer,intent(in) :: nlayer            ! Number of atmospheric layers.
215      integer,intent(in) :: nq                ! Number of tracers.
216
217      logical,intent(in) :: firstcall ! Signals first call to physics.
218      logical,intent(in) :: lastcall  ! Signals last call to physics.
219
220      real,intent(in) :: pday                  ! Number of elapsed sols since reference Ls=0.
221      real,intent(in) :: ptime                 ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
222      real,intent(in) :: ptimestep             ! Physics timestep (s).
223      real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
224      real,intent(in) :: pplay(ngrid,nlayer)   ! Mid-layer pressure (Pa).
225      real,intent(in) :: pphi(ngrid,nlayer)    ! Geopotential at mid-layer (m2s-2).
226      real,intent(in) :: pu(ngrid,nlayer)      ! Zonal wind component (m/s).
227      real,intent(in) :: pv(ngrid,nlayer)      ! Meridional wind component (m/s).
228      real,intent(in) :: pt(ngrid,nlayer)      ! Temperature (K).
229      real,intent(in) :: pq(ngrid,nlayer,nq)   ! Tracers (kg/kg_of_air).
230      real,intent(in) :: flxw(ngrid,nlayer)    ! Vertical mass flux (ks/s) at lower boundary of layer
231
232!   OUTPUTS:
233!   --------
234
235!     Physical tendencies :
236
237      real,intent(out) :: pdu(ngrid,nlayer)    ! Zonal wind tendencies (m/s/s).
238      real,intent(out) :: pdv(ngrid,nlayer)    ! Meridional wind tendencies (m/s/s).
239      real,intent(out) :: pdt(ngrid,nlayer)    ! Temperature tendencies (K/s).
240      real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
241      real,intent(out) :: pdpsrf(ngrid)        ! Surface pressure tendency (Pa/s).
242
243! Local saved variables:
244! ----------------------
245      integer,save :: day_ini                                      ! Initial date of the run (sol since Ls=0).
246      integer,save :: icount                                       ! Counter of calls to physiq during the run.
247!$OMP THREADPRIVATE(day_ini,icount)
248
249      !Pluto specific
250      REAL,save  :: acond,bcond
251      REAL,save  :: tcond1p4Pa
252      DATA tcond1p4Pa/38/
253
254! Local variables :
255! -----------------
256      !     Tendencies for the paleoclimate mode
257      REAL qsurfyear(ngrid,nq)  ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
258      REAL phisfinew(ngrid)       ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
259      REAL qsurfpal(ngrid,nq)           ! qsurf after a paleoclimate step : for physdem1 and restartfi
260      REAL phisfipal(ngrid)               ! geopotential after a paleoclimate step : for physdem1 and restartfi
261      REAL oblipal                   ! change of obliquity
262      REAL peri_daypal               ! new periday
263      REAL eccpal                    ! change of eccentricity
264      REAL tpalnew                   ! change of time
265      REAL adjustnew                 ! change in N2 ice albedo
266      REAL pdaypal                   ! new pday = day_ini + step
267      REAL zdt_tot                   ! time range corresponding to the flux of qsurfyear
268      REAL massacc(nq)             ! accumulated mass
269      REAL masslost(nq)            ! accumulated mass
270
271      REAL globave                   ! globalaverage 2D ps
272      REAL globaveice(nq)          ! globalaverage 2D ice
273      REAL globavenewice(nq)          ! globalaverage 2D ice
274      INTEGER lecttsoil     ! lecture of tsoil from proftsoil
275      REAL qsurf1(ngrid,nq)      ! saving qsurf to calculate flux over long timescales kg.m-2
276      REAL flusurf(ngrid,nq)     ! flux cond/sub kg.m-2.s-1
277      REAL flusurfold(ngrid,nq)  ! old flux cond/sub kg.m-2.s-1
278      REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
279
280      REAL,SAVE :: ptime0    ! store the first time
281      REAL dstep
282      REAL,SAVE :: glastep=20   ! step in pluto day to spread glacier
283
284!     Aerosol (dust or ice) extinction optical depth  at reference wavelength
285!     for the "naerkind" optically active aerosols:
286
287      real,save,allocatable :: dtau_aer(:,:,:) ! Aerosols
288!$OMP THREADPRIVATE(dtau_aer)
289      real zh(ngrid,nlayer)               ! Potential temperature (K).
290      real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
291      real omega(ngrid,nlayer)            ! omega velocity (Pa/s, >0 when downward)
292
293      integer i,l,ig,ierr,iq,nw,isoil,iesp
294
295      real zls                       ! Solar longitude (radians).
296      real zlss                      ! Sub solar point longitude (radians).
297      real zday                      ! Date (time since Ls=0, calculated in sols).
298      REAL,save :: saveday           ! saved date
299      REAL,save :: savedeclin        ! saved declin
300      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers.
301      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
302
303! TENDENCIES due to various processes :
304
305      ! For Surface Temperature : (K/s)
306      real zdtsurf(ngrid)     ! Cumulated tendencies.
307      real zdtsurfmr(ngrid)   ! Mass_redistribution routine.
308      real zdtsurfc(ngrid)    ! Condense_n2 routine.
309      real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
310
311      ! For Atmospheric Temperatures : (K/s)
312      real dtlscale(ngrid,nlayer)                             ! Largescale routine.
313      real zdtc(ngrid,nlayer)                                 ! Condense_n2 routine.
314      real zdtdif(ngrid,nlayer)                               ! Turbdiff/vdifc routines.
315      real zdtmr(ngrid,nlayer)                                ! Mass_redistribution routine.
316      real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
317      real zdtchim(ngrid,nlayer)                              ! Calchim routine.
318
319      ! For Surface Tracers : (kg/m2/s)
320      real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
321      !real zdqsurfc(ngrid)                  ! Condense_n2 routine.
322      REAL zdqsc(ngrid,nq)                  ! Condense_n2 routine.
323      real zdqsdif(ngrid,nq)                ! Turbdiff/vdifc routines.
324      real zdqssed(ngrid,nq)                ! Callsedim routine.
325      real zdqsurfmr(ngrid,nq)              ! Mass_redistribution routine.
326
327      ! For Tracers : (kg/kg_of_air/s)
328      real zdqc(ngrid,nlayer,nq)      ! Condense_n2 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.
332      real zdqsed(ngrid,nlayer,nq)    ! Callsedim routine.
333      real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
334      REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine
335      REAL,allocatable,save :: zdqschim(:,:)  ! Calchim_asis routine
336!$OMP THREADPRIVATE(zdqchim,zdqschim)
337
338      !! PLUTO variables
339      REAL zdqch4cloud(ngrid,nlayer,nq)
340      REAL zdqsch4cloud(ngrid,nq)
341      REAL zdtch4cloud(ngrid,nlayer)
342      REAL zdqcocloud(ngrid,nlayer,nq)
343      REAL zdqscocloud(ngrid,nq)
344      REAL zdtcocloud(ngrid,nlayer)
345      REAL rice_ch4(ngrid,nlayer)    ! Methane ice geometric mean radius (m)
346      REAL rice_co(ngrid,nlayer)     ! CO ice geometric mean radius (m)
347
348      REAL zdqsch4fast(ngrid)    ! used only for fast model nogcm
349      REAL zdqch4fast(ngrid)    ! used only for fast model nogcm
350      REAL zdqscofast(ngrid)    ! used only for fast model nogcm
351      REAL zdqcofast(ngrid)    ! used only for fast model nogcm
352      REAL zdqflow(ngrid,nq)
353
354      REAL zdtconduc(ngrid,nlayer) ! (K/s)
355      REAL zdumolvis(ngrid,nlayer)
356      REAL zdvmolvis(ngrid,nlayer)
357      real zdqmoldiff(ngrid,nlayer,nq)
358
359      ! Haze relatated tendancies
360      REAL zdqhaze(ngrid,nlayer,nq)
361      REAL zdqprodhaze(ngrid,nq)
362      REAL zdqsprodhaze(ngrid)
363      REAL zdqhaze_col(ngrid)
364      REAL zdqphot_prec(ngrid,nlayer)
365      REAL zdqphot_ch4(ngrid,nlayer)
366      REAL zdqconv_prec(ngrid,nlayer)
367      REAL zdq_source(ngrid,nlayer,nq)
368      ! Fast Haze relatated tendancies
369      REAL fluxbot(ngrid)
370      REAL gradflux(ngrid)
371      REAL fluxlym_sol_bot(ngrid)      ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
372      REAL fluxlym_ipm_bot(ngrid)      ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
373      REAL flym_sol(ngrid)      ! Incident Solar flux Lyman alpha ph.m-2.s-1
374      REAL flym_ipm(ngrid)      ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
375      REAL zfluxuv                     ! Lyman alpha flux at 1AU
376
377      REAL array_zero1(ngrid)
378      REAL array_zero2(ngrid,nlayer)
379
380      ! For Winds : (m/s/s)
381      real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer)       ! Convadj routine.
382      real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer)         ! Mass_redistribution routine.
383      real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer)       ! Turbdiff/vdifc routines.
384      real zdhdif(ngrid,nlayer)                             ! Turbdiff/vdifc routines.
385      real zdhadj(ngrid,nlayer)                             ! Convadj routine.
386      REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)            ! condense_n2 routine.
387
388      ! For Pressure and Mass :
389      real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
390      real zdmassmr_col(ngrid)    ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
391      real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
392
393      ! Local variables for MICROPHYSICS:
394      ! ---------------------------------
395      real gzlat(ngrid,nlayer)           ! Altitude-Latitude-dependent gravity (this should be stored elsewhere...).
396      real pdqmufi(ngrid,nlayer,nq)      ! Microphysical tendency (X/kg_of_air/s).
397      real pdqmufi_prod(ngrid,nlayer,nq) ! Aerosols production tendency (kg/kg_of_air/s).
398      real int2ext(ngrid,nlayer)         ! Intensive to extensive factor (kg_air/m3: X/kg_air --> X/m3).
399
400! Local variables for LOCAL CALCULATIONS:
401! ---------------------------------------
402      real zflubid(ngrid)
403      real zplanck(ngrid),zpopsk(ngrid,nlayer)
404      REAL zdum1(ngrid,nlayer)
405      REAL zdum2(ngrid,nlayer)
406      real ztim1,ztim2,ztim3, z1,z2
407      real ztime_fin
408      real zdh(ngrid,nlayer)
409      real gmplanet
410      real taux(ngrid),tauy(ngrid)
411
412! local variables for DIAGNOSTICS : (diagfi & stat)
413! -------------------------------------------------
414      real ps(ngrid)                                     ! Surface Pressure.
415      real zt(ngrid,nlayer)                              ! Atmospheric Temperature.
416      real zu(ngrid,nlayer),zv(ngrid,nlayer)             ! Zonal and Meridional Winds.
417      real zq(ngrid,nlayer,nq)                           ! Atmospheric Tracers.
418      real zdtadj(ngrid,nlayer)                          ! Convadj Diagnostic.
419      real zdtdyn(ngrid,nlayer)                          ! Dynamical Heating (K/s).
420      real zdudyn(ngrid,nlayer)                          ! Dynamical Zonal Wind tendency (m.s-2).
421
422      real reff(ngrid,nlayer)                       ! Effective dust radius (used if doubleq=T).
423      real vmr(ngrid,nlayer)                        ! volume mixing ratio
424      real time_phys
425
426      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
427
428      real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
429
430      !     Pluto adding variables
431      real vmr_ch4(ngrid)  ! vmr ch4
432      real vmr_co(ngrid)  ! vmr co
433      real rho(ngrid,nlayer) ! density
434      real zrho_ch4(ngrid,nlayer) ! density methane kg.m-3
435      real zrho_co(ngrid,nlayer) ! density CO kg.m-3
436      real zrho_haze(ngrid,nlayer) ! density haze kg.m-3
437      real zdqrho_photprec(ngrid,nlayer) !photolysis rate kg.m-3.s-1
438      real zq1temp_ch4(ngrid) !
439      real qsat_ch4(ngrid) !
440      real qsat_ch4_l1(ngrid) !
441!      CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
442      real profmmr(ngrid,nlayer) ! fixed profile of haze if haze_proffix
443      real sensiblehf1(ngrid) ! sensible heat flux
444      real sensiblehf2(ngrid) ! sensible heat flux
445
446!     included by RW for H2O Manabe scheme
447      real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
448      real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
449
450!     to test energy conservation (RW+JL)
451      real mass(ngrid,nlayer),massarea(ngrid,nlayer)
452      real dEtot, dEtots, AtmToSurf_TurbFlux
453      real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
454!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
455
456!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
457      real dtmoist_max,dtmoist_min
458      real dItot, dItot_tmp, dVtot, dVtot_tmp
459      real dWtot, dWtot_tmp, dWtots, dWtots_tmp
460
461      real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
462
463      real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
464
465!  Non-oro GW tendencies
466      REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
467      REAL d_t_hin(ngrid,nlayer)
468!  Diagnostics 2D of gw_nonoro
469      REAL zustrhi(ngrid), zvstrhi(ngrid)
470
471      real :: tsurf2(ngrid)
472!!      real :: flux_o(ngrid),flux_g(ngrid)
473      real :: flux_g(ngrid)
474      real :: flux_sens_lat(ngrid)
475      real :: qsurfint(ngrid,nq)
476
477      ! local variables for skin depth check
478      real :: therm_inertia(ngrid,nsoilmx)
479    !   real :: tidat_out(ngrid,nsoilmx)
480      real :: inertia_min,inertia_max
481      real :: diurnal_skin ! diurnal skin depth (m)
482      real :: annual_skin ! anuual skin depth (m)
483
484      ! when no startfi file is asked for init
485      real alpha,lay1 ! coefficients for building layers
486      integer iloop
487
488      ! flags to trigger extra sanity checks
489      logical, save ::  check_physics_inputs=.false.
490      logical, save ::  check_physics_outputs=.false.
491!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
492
493      ! Misc
494      character*2 :: str2
495      character(len=10) :: tmp1
496      character(len=10) :: tmp2
497!==================================================================================================
498
499! -----------------
500! I. INITIALISATION
501! -----------------
502
503! --------------------------------
504! I.1   First Call Initialisation.
505! --------------------------------
506      if (firstcall) then
507         call getin_p("check_physics_inputs", check_physics_inputs)
508         call getin_p("check_physics_outputs", check_physics_outputs)
509
510         ! Allocate saved arrays (except for 1D model, where this has already
511         ! been done)
512         if (ngrid>1) call phys_state_var_init(nq)
513
514!        Variables set to 0
515!        ~~~~~~~~~~~~~~~~~~
516         dtrad(:,:) = 0.0
517         fluxrad(:) = 0.0
518         tau_col(:) = 0.0
519         zdtsw(:,:) = 0.0
520         zdtlw(:,:) = 0.0
521
522!        Initialize tracer names, indexes and properties.
523!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
524         IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
525         if (tracer) then
526            call initracer(ngrid,nq)
527            ! if(photochem) then !AF24: removed
528         endif
529!        Initialize aerosol indexes.
530!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~
531         ! call iniaerosol
532         ! allocate related local arrays
533         ! (need be allocated instead of automatic because of "naerkind")
534         allocate(dtau_aer(ngrid,nlayer,naerkind))
535
536#ifdef CPP_XIOS
537         ! Initialize XIOS context
538         write(*,*) "physiq: call wxios_context_init"
539         CALL wxios_context_init
540#endif
541
542!        Read 'startfi.nc' file.
543!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
544         call phyetat0(startphy_file,                                 &
545                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
546                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
547
548         if (.not.startphy_file) then
549           ! starting without startfi.nc and with callsoil
550           ! is not yet possible as soildepth default is not defined
551           if (callsoil) then
552              ! default mlayer distribution, following a power law:
553              !  mlayer(k)=lay1*alpha**(k-1/2)
554              lay1=2.e-4
555              alpha=2
556              do iloop=0,nsoilmx-1
557                 mlayer(iloop)=lay1*(alpha**(iloop-0.5))
558              enddo
559              lay1=sqrt(mlayer(0)*mlayer(1))
560              alpha=mlayer(1)/mlayer(0)
561              do iloop=1,nsoilmx
562                 layer(iloop)=lay1*(alpha**(iloop-1))
563              enddo
564           endif
565           ! additionnal "academic" initialization of physics
566           if (is_master) write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
567           tsurf(:)=pt(:,1)
568           if (is_master) write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!"
569           do isoil=1,nsoilmx
570             tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
571           enddo
572           if (is_master) write(*,*) "Physiq: initializing day_ini to pday !"
573           day_ini=pday
574         endif
575
576         if (pday.ne.day_ini) then
577             write(*,*) "ERROR in physiq.F90:"
578             write(*,*) "bad synchronization between physics and dynamics"
579             write(*,*) "dynamics day: ",pday
580             write(*,*) "physics day:  ",day_ini
581             stop
582         endif
583
584         write (*,*) 'In physiq day_ini =', day_ini
585
586!        Initialize albedo calculation.
587!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
588         albedo(:,:)=0.0
589         albedo_bareground(:)=0.0
590         albedo_snow_SPECTV(:)=0.0
591         albedo_n2_ice_SPECTV(:)=0.0
592
593         ptime0=ptime
594         write (*,*) 'In physiq ptime0 =', ptime
595
596         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
597
598!        Initialize orbital calculation.
599!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
600         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
601
602         savedeclin=0.
603         saveday=pday
604         adjust=0. ! albedo adjustment for convergeps
605
606!        Initialize soil.
607!        ~~~~~~~~~~~~~~~~
608         if (callsoil) then
609            call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
610                      ptimestep,tsurf,tsoil,capcal,fluxgrd)
611         else ! else of 'callsoil'.
612            print*,'WARNING! Thermal conduction in the soil turned off'
613            capcal(:)=1.e6
614            fluxgrd(:)=intheat
615            print*,'Flux from ground = ',intheat,' W m^-2'
616         endif ! end of 'callsoil'.
617
618         icount=1
619
620!        Initialize variable for dynamical heating and zonal wind tendency diagnostic
621!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
622         ztprevious(:,:)=pt(:,:)
623         zuprevious(:,:)=pu(:,:)
624
625         if(meanOLR)then
626            call system('rm -f rad_bal.out') ! to record global radiative balance.
627            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
628            call system('rm -f h2o_bal.out') ! to record global hydrological balance.
629         endif
630
631         if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
632            call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
633                         ptimestep,pday+nday,time_phys,cell_area,          &
634                         albedo_bareground,zmea,zstd,zsig,zgam,zthe)
635         endif
636
637!        Initialize correlated-k.
638!        ~~~~~~~~~~~~~~~~~~~~~~~~
639         if (corrk) then
640            ! We initialise the spectral grid here instead of
641            ! at firstcall of callcorrk so we can output XspecIR, XspecVI
642            ! when using Dynamico
643            if (is_master) print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
644            call getin_p("corrkdir",corrkdir)
645            if (is_master) print*,"corrkdir = ", corrkdir
646            write (tmp1, '(i4)') L_NSPECTI
647            write (tmp2, '(i4)') L_NSPECTV
648            banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)))
649            banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)))
650            call setspi !Basic infrared properties.
651            call setspv ! Basic visible properties.
652            call sugas_corrk       ! Set up gaseous absorption properties.
653            if (optichaze) then
654               call suaer_corrk       ! Set up aerosol optical properties.
655            endif
656         endif
657
658!        Initialize microphysics.
659!        ~~~~~~~~~~~~~~~~~~~~~~~~
660         IF (callmufi) THEN
661            ! Initialize microphysics arrays.
662            call inimufi(ptimestep)
663         ENDIF ! end callmufi
664
665!!         call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
666         ! XIOS outputs
667#ifdef CPP_XIOS
668
669         if (is_master) write(*,*) "physiq: call initialize_xios_output"
670         call initialize_xios_output(pday,ptime,ptimestep,daysec,year_day, &
671                                     presnivs,pseudoalt,mlayer,WNOI,WNOV)
672#endif
673
674!!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
675
676         write(*,*) "physiq: end of firstcall"
677      endif ! end of 'firstcall'
678
679!!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
680!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
681
682      if (check_physics_inputs) then
683         !check the validity of input fields coming from the dynamics
684         call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
685      endif
686
687!      call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
688!      call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
689
690! ------------------------------------------------------
691! I.2   Initializations done at every physical timestep:
692! ------------------------------------------------------
693
694#ifdef CPP_XIOS
695      ! update XIOS time/calendar
696      call update_xios_timestep
697#endif
698
699      ! Initialize various variables
700      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
701
702      pdt(1:ngrid,1:nlayer) = 0.0
703      zdtsurf(1:ngrid)      = 0.0
704      pdq(1:ngrid,1:nlayer,1:nq) = 0.0
705      dqsurf(1:ngrid,1:nq)= 0.0
706      pdu(1:ngrid,1:nlayer) = 0.0
707      pdv(1:ngrid,1:nlayer) = 0.0
708      pdpsrf(1:ngrid)       = 0.0
709      zflubid(1:ngrid)      = 0.0
710      flux_sens_lat(1:ngrid) = 0.0
711      taux(1:ngrid) = 0.0
712      tauy(1:ngrid) = 0.0
713
714      if (conservn2) then
715         write(*,*) 'conservn2 iniloop'
716         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
717      endif
718
719      zday=pday+ptime ! Compute time, in sols (and fraction thereof).
720
721      ! Compute Stellar Longitude (Ls), and orbital parameters.
722      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
723      if (season) then
724         call stellarlong(zday,zls)
725      else
726         call stellarlong(noseason_day,zls)
727      end if
728
729!     Get Lyman alpha flux at specific Ls
730      if (callmufi.or.haze) then
731         call lymalpha(zls,zfluxuv)
732         print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
733      end if
734
735      IF (triton) then
736         CALL orbitetriton(zls,zday,dist_star,declin)
737      ELSE
738         call orbite(zls,dist_star,declin,right_ascen)
739      ENDIF
740
741      if (diurnal) then
742              zlss=-2.*pi*(zday-.5)
743      else if(diurnal .eqv. .false.) then
744              zlss=9999.
745      endif
746
747      glat(:) = g !AF24: removed oblateness
748
749      ! Compute geopotential between layers.
750      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
751      zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
752      do l=1,nlayer
753         zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
754      enddo
755
756      zzlev(1:ngrid,1)=0.
757
758      do l=2,nlayer
759         do ig=1,ngrid
760            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
761            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
762            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
763         enddo
764      enddo
765
766      !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
767
768      zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
769
770      ! Compute potential temperature
771      ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
772      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
773      do l=1,nlayer
774         do ig=1,ngrid
775            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
776            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
777            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
778            massarea(ig,l)=mass(ig,l)*cell_area(ig)
779         enddo
780      enddo
781
782     ! Compute vertical velocity (m/s) from vertical mass flux
783     ! w = F / (rho*area) and rho = P/(r*T)
784     ! But first linearly interpolate mass flux to mid-layers
785      if (.not.fast) then
786       do l=1,nlayer-1
787         pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
788       enddo
789       pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
790       do l=1,nlayer
791         pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /  &
792                       (pplay(1:ngrid,l)*cell_area(1:ngrid))
793       enddo
794       ! omega in Pa/s
795       do l=1,nlayer-1
796         omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
797       enddo
798       omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
799       do l=1,nlayer
800         omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
801       enddo
802      endif
803
804      if (conservn2) then
805         write(*,*) 'conservn2 thermo'
806         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
807      endif
808
809      !  Compute variations of g with latitude (to do).
810      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
811      gzlat(:,:) = g
812
813      ! Initialize microphysical diagnostics.
814      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
815      IF (callmufi) THEN
816         ! Initialize intensive to extensive factor (kg_air/m3: X/kg_air --> X/m3).
817         int2ext(:,:) = (pplev(:,1:nlayer)-pplev(:,2:nlayer+1)) / gzlat(:,1:nlayer) / (zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
818
819         ! Initialize microphysics diagnostics arrays.
820         call inimufi_diag(ngrid,nlayer,nq,pq,int2ext)
821      ENDIF ! end callmufi
822
823! --------------------------------------------------------
824!    II.1 Thermosphere
825! --------------------------------------------------------
826
827! ajout de la conduction depuis la thermosphere
828      IF (callconduct) THEN
829
830          call conduction (ngrid,nlayer,ptimestep, &
831                      pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
832          DO l=1,nlayer
833             DO ig=1,ngrid
834                pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
835             ENDDO
836          ENDDO
837
838      ENDIF
839
840! ajout de la viscosite moleculaire
841      IF (callmolvis) THEN
842          call molvis(ngrid,nlayer,ptimestep,   &
843                      pplay,pt,zzlay,zzlev,  &
844                      zdtconduc,pu,tsurf,zdumolvis)
845          call molvis(ngrid,nlayer,ptimestep,   &
846                      pplay,pt,zzlay,zzlev,  &
847                      zdtconduc,pv,tsurf,zdvmolvis)
848
849          DO l=1,nlayer
850             DO ig=1,ngrid
851             ! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
852                pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
853                pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
854             ENDDO
855          ENDDO
856      ENDIF
857
858      IF (callmoldiff) THEN
859           call moldiff_red(ngrid,nlayer,nq, &
860                        pplay,pplev,pt,pdt,pq,pdq,ptimestep,   &
861                        zzlay,zdtconduc,zdqmoldiff)
862
863           DO l=1,nlayer
864              DO ig=1,ngrid
865                 DO iq=1, nq
866                  pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
867                 ENDDO
868              ENDDO
869           ENDDO
870      ENDIF
871
872      if (conservn2) then
873         write(*,*) 'conservn2 thermosphere'
874         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
875      endif
876
877
878!---------------------------------
879! II.2 Compute radiative tendencies
880!---------------------------------
881!     Saving qsurf to compute paleo flux condensation/sublimation
882      DO iq=1, nq
883         DO ig=1,ngrid
884             IF (qsurf(ig,iq).lt.0.) then
885                qsurf(ig,iq)=0.
886             ENDIF
887             qsurf1(ig,iq)=qsurf(ig,iq)
888         ENDDO
889      ENDDO
890
891
892      ! Compute local stellar zenith angles
893      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
894      fract = 0
895      if (diurnal) then
896         ztim1=SIN(declin)
897         ztim2=COS(declin)*COS(2.*pi*(zday-.5))
898         ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
899
900         call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
901                        ztim1,ztim2,ztim3,mu0,fract)
902      else if(diurnal .eqv. .false.) then
903
904         call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
905         ! WARNING: this function appears not to work in 1D
906
907         if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
908            mu0 = cos(pi*szangle/180.0)
909            fract= 1/(4*mu0) ! AF24: from pluto.old
910         endif
911
912      endif
913
914
915!     Pluto albedo /IT changes depending on surface ices (only in 3D)
916!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
917      if (ngrid.ne.1) then
918
919         !! Specific to change albedo of N2 so that Psurf
920         !! converges toward 1.4 Pa in "1989" seasons for Triton
921         !! converges toward 1.1 Pa in "2015" seasons for Pluto
922         if (convergeps) then
923            if (triton) then
924               ! 1989 declination
925               if (declin*180./pi.gt.-46..and.declin*180./pi.lt.-45.   &
926            .and.zday.gt.saveday+1000.   &
927            .and.declin.lt.savedeclin) then
928               call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
929               if (globave.gt.1.5) then
930                     adjust=adjust+0.005
931               else if (globave.lt.1.3) then
932                     adjust=adjust-0.005
933               endif
934               saveday=zday
935               endif
936            else
937               ! Pluto : 2015 declination current epoch
938               if (declin*180./pi.gt.50.and.declin*180./pi.lt.51.  &
939            .and.zday.gt.saveday+10000.  &
940            .and.declin.gt.savedeclin) then
941               call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
942               if (globave.gt.1.2) then
943                     adjust=adjust+0.005
944               else if (globave.lt.1.) then
945                     adjust=adjust-0.005
946               endif
947               saveday=zday
948               endif
949            endif
950         end if
951      end if ! if ngrid ne 1
952
953      call surfprop(ngrid,nq,fract,qsurf,tsurf,        &
954      capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
955      albedo,emis,therm_inertia)
956      ! do i=2,ngrid
957      !    albedo(i,:) = albedo(1,:)
958      ! enddo
959      ! AF24: TODO check albedo has been initialized here
960
961      if (firstcall.and.callsoil) then
962         ! AF24 Originally in soil.F, but therm_inertia is modified by surfprop
963         ! Additional checks: is the vertical discretization sufficient
964         ! to resolve diurnal and annual waves?
965         do ig=1,ngrid
966            ! extreme inertia for this column
967            inertia_min=minval(therm_inertia(ig,:))
968            inertia_max=maxval(therm_inertia(ig,:))
969            ! diurnal and annual skin depth
970            diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
971            annual_skin=(inertia_max/volcapa)*sqrt(year_day*daysec/pi)
972            if (0.5*diurnal_skin<layer(1)) then
973            ! one should have the fist layer be at least half of diurnal skin depth
974            write(*,*) "soil Error: grid point ig=",ig
975            write(*,*) "            longitude=",longitude(ig)*(180./pi)
976            write(*,*) "             latitude=",latitude(ig)*(180./pi)
977            write(*,*) "  first soil layer depth ",layer(1)
978            write(*,*) "  not small enough for a diurnal skin depth of ", &
979                        diurnal_skin
980            write(*,*) " change soil layer distribution (comsoil_h.F90)"
981            call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
982            endif
983            if (2.*annual_skin>layer(nsoilmx)) then
984            ! one should have the full soil be at least twice the diurnal skin depth
985            write(*,*) "soil Error: grid point ig=",ig
986            write(*,*) "            longitude=",longitude(ig)*(180./pi)
987            write(*,*) "             latitude=",latitude(ig)*(180./pi)
988            write(*,*) "  total soil layer depth ",layer(nsoilmx)
989            write(*,*) "  not large enough for an annual skin depth of ", &
990                        annual_skin
991            write(*,*) " change soil layer distribution (comsoil_h.F90)"
992            call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
993            endif
994         enddo ! of do ig=1,ngrid
995
996      end if ! callsoil
997
998      if (callrad) then
999         if( mod(icount-1,iradia).eq.0.or.lastcall) then
1000
1001            ! Eclipse incoming sunlight !AF24: removed
1002
1003!!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1004!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1005
1006
1007            if (corrk) then
1008
1009! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1010! II.a Call correlated-k radiative transfer scheme
1011! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1012               if(kastprof)then
1013                  print*,'kastprof should not = true here'
1014                  call abort
1015               endif
1016
1017               ! standard callcorrk
1018               if (oldplutocorrk) then
1019                  call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf,          &
1020                               albedo(:,1),emis,mu0,pplev,pplay,pt,              &
1021                               zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,       &
1022                               zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,   &
1023                               fluxtop_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
1024                               firstcall,lastcall)
1025                  albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1026                  fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+       &
1027                               fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid,1))
1028                  fluxabs_sw(1:ngrid)=fluxtop_dn(1:ngrid)-fluxtop_sw(1:ngrid)
1029               else
1030                muvar(1:ngrid,1:nlayer+1)=mugaz
1031                call callcorrk(ngrid,nlayer,pq,nq,qsurf,  &
1032                              albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
1033                              zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,muvar,   &
1034                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
1035                              fluxsurfabs_sw,fluxtop_lw,                          &
1036                              fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu,         &
1037                              int_dtaui,int_dtauv,                                &
1038                              tau_col,firstcall,lastcall)
1039                  ! Radiative flux from the sky absorbed by the surface (W.m-2).
1040                  GSR=0.0
1041                  fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+        &
1042                              fluxsurfabs_sw(1:ngrid)
1043               endif ! oldplutocorrk
1044                !GG (feb2021): Option to "artificially" decrease the raditive time scale in
1045                !the deep atmosphere  press > 0.1 bar. Suggested by MT.
1046                !! COEFF_RAD to be "tuned" to facilitate convergence of tendency
1047
1048                !coeff_rad=0.   ! 0 values, it doesn't accelerate the convergence
1049                !coeff_rad=0.5
1050                !coeff_rad=1.
1051                !do l=1, nlayer
1052                !  do ig=1,ngrid
1053                !    if(pplay(ig,l).ge.1.d4) then
1054                !      zdtsw(ig,l)=zdtsw(ig,l)*(pplay(ig,l)/1.d4)**coeff_rad
1055                !      zdtlw(ig,l)=zdtlw(ig,l)*(pplay(ig,l)/1.d4)**coeff_rad
1056                !    endif
1057                !  enddo
1058                !enddo
1059
1060                ! AF24: removed CLFvarying Option
1061
1062
1063                            !if(noradsurf)then ! no lower surface; SW flux just disappears
1064                            !   GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
1065                            !   fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
1066                            !   print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
1067                            !endif
1068
1069               ! Net atmospheric radiative heating rate (K.s-1)
1070               dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
1071
1072               ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
1073               if (firstcall .and. albedo_spectral_mode) then
1074                  call spectral_albedo_calc(albedo_snow_SPECTV)
1075               endif
1076
1077            else
1078! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1079! II.b Atmosphere has no radiative effect
1080! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1081               fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
1082               if(ngrid.eq.1)then ! / by 4 globally in 1D case...
1083                  fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
1084               endif
1085               fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
1086               print*,'------------WARNING---WARNING------------' ! by MT2015.
1087               print*,'You are in corrk=false mode, '
1088               print*,'and the surface albedo is taken equal to the first visible spectral value'
1089
1090               albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
1091               fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
1092               fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid)
1093               fluxrad_sky(1:ngrid)    = fluxsurfabs_sw(1:ngrid)
1094               fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
1095
1096               dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
1097
1098            endif ! end of corrk
1099
1100         endif ! of if(mod(icount-1,iradia).eq.0)
1101
1102
1103         ! Transformation of the radiative tendencies
1104         ! ------------------------------------------
1105         zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
1106         zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid)
1107         fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
1108         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
1109
1110         ! Test of energy conservation
1111         !----------------------------
1112         if(enertest)then
1113            call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
1114            call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW)
1115            !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
1116            call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
1117            call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*cell_area(:)/totarea_planet,dEtotsLW)
1118            dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
1119            dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
1120            if (is_master) then
1121               print*,'---------------------------------------------------------------'
1122               print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
1123               print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
1124               print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
1125               print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
1126               print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
1127               print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
1128            endif
1129         endif ! end of 'enertest'
1130
1131      endif ! of if (callrad)
1132
1133!!      call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1134!!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1135
1136   if (conservn2) then
1137      write(*,*) 'conservn2 radiat'
1138      call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
1139   endif
1140
1141!  --------------------------------------------
1142!  III. Vertical diffusion (turbulent mixing) :
1143!  --------------------------------------------
1144
1145      if (calldifv) then
1146
1147         zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
1148
1149         if (oldplutovdifc) then
1150            zdum1(:,:) = 0.
1151            zdum2(:,:) = 0.
1152            zdh(:,:)=pdt(:,:)/zpopsk(:,:)
1153
1154            ! Calling vdif (Martian version WITH N2 condensation)
1155            CALL vdifc_pluto(ngrid,nlayer,nq,zpopsk,       &
1156                    ptimestep,capcal,lwrite,         &
1157                    pplay,pplev,zzlay,zzlev,z0,      &
1158                    pu,pv,zh,pq,pt,tsurf,emis,qsurf, &
1159                    zdum1,zdum2,zdh,pdq,pdt,zflubid, &
1160                    zdudif,zdvdif,zdhdif,zdtsdif,q2, &
1161                    zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
1162
1163            zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1164
1165            bcond=1./tcond1p4Pa
1166            acond=r/lw_n2
1167
1168         ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
1169         else if (UseTurbDiff) then
1170
1171            call turbdiff(ngrid,nlayer,nq,                  &
1172                          ptimestep,capcal,                      &
1173                          pplay,pplev,zzlay,zzlev,z0,            &
1174                          pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
1175                          pdt,pdq,zflubid,                       &
1176                          zdudif,zdvdif,zdtdif,zdtsdif,          &
1177                          sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
1178                          taux,tauy)
1179
1180         else ! if .not. (oldplutovdifc) .and. (UseTurbDiff)
1181
1182            zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1183
1184            call vdifc(ngrid,nlayer,nq,zpopsk,           &
1185                       ptimestep,capcal,lwrite,               &
1186                       pplay,pplev,zzlay,zzlev,z0,            &
1187                       pu,pv,zh,pq,tsurf,emis,qsurf,          &
1188                       zdh,pdq,zflubid,                       &
1189                       zdudif,zdvdif,zdhdif,zdtsdif,          &
1190                       sensibFlux,q2,zdqdif,zdqsdif)
1191
1192            zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1193            zdqevap(1:ngrid,1:nlayer)=0.
1194
1195         end if !end of 'UseTurbDiff'
1196
1197         zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
1198
1199         !!! this is always done, except for turbulence-resolving simulations
1200         if (.not. turb_resolved) then
1201           pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
1202           pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
1203           pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
1204         endif
1205
1206         if (tracer) then
1207           pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
1208           dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
1209         end if ! of if (tracer)
1210
1211         ! test energy conservation
1212         !-------------------------
1213         if(enertest)then
1214
1215            dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
1216            do ig = 1, ngrid
1217               dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
1218               dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
1219            enddo
1220
1221            call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
1222            dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
1223            call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
1224            call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
1225
1226            if (is_master) then
1227
1228               if (UseTurbDiff) then
1229                         print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1230                         print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
1231                  print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1232               else
1233                         print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
1234                         print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
1235                         print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1236               end if
1237            endif ! end of 'is_master'
1238
1239         ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
1240         endif ! end of 'enertest'
1241
1242      else ! calldifv
1243
1244         ztim1=4.*sigma*ptimestep
1245         DO ig=1,ngrid
1246           ztim2=ztim1*emis(ig)*tsurf(ig)**3
1247           z1=capcal(ig)*tsurf(ig)+  &
1248           ztim2*tsurf(ig)+ (fluxrad(ig)+fluxgrd(ig))*ptimestep
1249           z2= capcal(ig)+ztim2
1250           zdtsurf(ig)=(z1/z2 - tsurf(ig))/ptimestep
1251
1252           ! for output:
1253           !dplanck(ig)=4.*stephan*ptimestep*emis(ig)*tsurf(ig)**3
1254         ENDDO
1255
1256         ! if(.not.newtonian)then
1257         !zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
1258
1259!        ------------------------------------------------------------------
1260!        Methane surface sublimation and condensation in fast model (nogcm)
1261!        ------------------------------------------------------------------
1262         if ((methane).and.(fast).and.condmetsurf) THEN
1263
1264            call ch4surf(ngrid,nlayer,nq,ptimestep,capcal, &
1265               tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
1266               zdqch4fast,zdqsch4fast)
1267
1268            dqsurf(1:ngrid,igcm_ch4_ice)= dqsurf(1:ngrid,igcm_ch4_ice) + &
1269                                         zdqsch4fast(1:ngrid)
1270            pdq(1:ngrid,1,igcm_ch4_gas)= pdq(1:ngrid,1,igcm_ch4_gas) + &
1271                                         zdqch4fast(1:ngrid)
1272            zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_ch4*zdqsch4fast(1:ngrid)/capcal(1:ngrid)
1273            end if
1274!        ------------------------------------------------------------------
1275!        CO surface sublimation and condensation in fast model (nogcm)
1276!        ------------------------------------------------------------------
1277         if ((carbox).and.(fast).and.condcosurf) THEN
1278
1279            call cosurf(ngrid,nlayer,nq,ptimestep, &
1280               tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
1281               zdqcofast,zdqscofast)
1282
1283            dqsurf(1:ngrid,igcm_co_ice)= dqsurf(1:ngrid,igcm_co_ice) + &
1284                                         zdqscofast(1:ngrid)
1285            pdq(1:ngrid,1,igcm_co_gas)= pdq(1:ngrid,1,igcm_co_gas) + &
1286                                        zdqcofast(1:ngrid)
1287            zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_co*zdqscofast(1:ngrid)/capcal(1:ngrid)
1288         end if
1289
1290
1291      endif ! end of 'calldifv'
1292
1293      if (conservn2) then
1294        write(*,*) 'conservn2 calldifv'
1295        call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ &
1296                                       dqsurf(:,1)*ptimestep)
1297      endif
1298      if (methane.and.conservch4) then
1299        write(*,*) 'conservch4 calldifv'
1300        if (fast) then
1301             call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
1302                   qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
1303                   ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' vdifc ')
1304        else
1305             call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
1306                   igcm_ch4_gas,igcm_ch4_ice, &
1307                   ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
1308        endif
1309      endif
1310
1311!-------------------
1312!   IV. Convection :
1313!-------------------
1314
1315! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1316! IV.a Dry convective adjustment :
1317! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1318
1319      if(calladj) then
1320
1321         zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
1322         zduadj(1:ngrid,1:nlayer)=0.0
1323         zdvadj(1:ngrid,1:nlayer)=0.0
1324         zdhadj(1:ngrid,1:nlayer)=0.0
1325
1326
1327         call convadj(ngrid,nlayer,nq,ptimestep,            &
1328                      pplay,pplev,zpopsk,                   &
1329                      pu,pv,zh,pq,                          &
1330                      pdu,pdv,zdh,pdq,                      &
1331                      zduadj,zdvadj,zdhadj,                 &
1332                      zdqadj)
1333
1334         pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
1335         pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
1336         pdt(1:ngrid,1:nlayer)    = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
1337         zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
1338
1339         if(tracer) then
1340            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
1341         end if
1342
1343         ! Test energy conservation
1344         if(enertest)then
1345            call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot)
1346            if (is_master) print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
1347         endif
1348
1349         ! ! Test water conservation !AF24: removed
1350
1351      endif ! end of 'calladj'
1352
1353!-----------------------------------------------
1354!   V. Nitrogen condensation-sublimation :
1355!-----------------------------------------------
1356
1357      if (n2cond) then
1358
1359         if (.not.tracer) then
1360            print*,'We need a N2 ice tracer to condense N2'
1361            call abort
1362         endif
1363
1364         zdqc(:,:,:)=0.
1365         zdqsc(:,:)=0.
1366         call condense_n2(ngrid,nlayer,nq,ptimestep,                    &
1367                           capcal,pplay,pplev,tsurf,pt,                 &
1368                           pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,       &
1369                           qsurf(1,igcm_n2),albedo,emis,                &
1370                           zdtc,zdtsurfc,pdpsrf,zduc,zdvc,              &
1371                           zdqc,zdqsc(1,igcm_n2))
1372
1373         pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
1374         pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
1375         pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
1376         zdtsurf(1:ngrid)      = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
1377
1378         pdq(1:ngrid,1:nlayer,1:nq)   = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
1379         dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
1380
1381!!         call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1382!!         call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
1383
1384         ! test energy conservation
1385         if(enertest)then
1386            call planetwide_sumval(cpp*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot)
1387            call planetwide_sumval(capcal(:)*zdtsurfc(:)*cell_area(:)/totarea_planet,dEtots)
1388            if (is_master) then
1389               print*,'In n2cloud atmospheric energy change   =',dEtot,' W m-2'
1390               print*,'In n2cloud surface energy change       =',dEtots,' W m-2'
1391            endif
1392         endif
1393
1394      endif  ! end of 'n2cond'
1395
1396      if (conservn2) then
1397       write(*,*) 'conservn2 n2cond'
1398       call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1399           pdpsrf(:)*ptimestep,qsurf(:,1)+dqsurf(:,1)*ptimestep)
1400      endif
1401      if (methane.and.conservch4) then
1402        write(*,*) 'conservch4 n2cond'
1403        if (fast) then
1404             call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
1405                   qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
1406                   ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' n2cond')
1407        else
1408             call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
1409                   igcm_ch4_gas,igcm_ch4_ice, &
1410                   ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond')
1411        endif
1412      endif
1413
1414!---------------------------------------------
1415!   VI. Specific parameterizations for tracers
1416!---------------------------------------------
1417
1418      if (tracer) then
1419
1420!      ---------------------------------------
1421!      Methane ice condensation in the atmosphere
1422!      ----------------------------------------
1423         rice_ch4(:,:)=0 ! initialization needed for callsedim
1424         zdqch4cloud(:,:,:)=0.
1425         if ((methane).and.(metcloud).and.(.not.fast)) THEN
1426            call ch4cloud(ngrid,nlayer,ptimestep,  &
1427                      pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,  &
1428                      pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud,   &
1429                      nq,rice_ch4)
1430
1431            DO l=1,nlayer
1432               DO ig=1,ngrid
1433                  pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+  &
1434                                         zdqch4cloud(ig,l,igcm_ch4_gas)
1435                  pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+  &
1436                                         zdqch4cloud(ig,l,igcm_ch4_ice)
1437               ENDDO
1438            ENDDO
1439
1440            ! Increment methane ice surface tracer tendency
1441            DO ig=1,ngrid
1442               dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+   &
1443                                     zdqsch4cloud(ig,igcm_ch4_ice)
1444            ENDDO
1445
1446            ! update temperature tendancy
1447            DO ig=1,ngrid
1448               DO l=1,nlayer
1449               pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
1450               ENDDO
1451            ENDDO
1452         end if
1453
1454!      ---------------------------------------
1455!      CO ice condensation in the atmosphere
1456!      ----------------------------------------
1457         zdqcocloud(:,:,:)=0.
1458         IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
1459            call cocloud(ngrid,nlayer,ptimestep,   &
1460                      pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,  &
1461                      pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud,   &
1462                      nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))
1463
1464            DO l=1,nlayer
1465               DO ig=1,ngrid
1466                  pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+ &
1467                                         zdqcocloud(ig,l,igcm_co_gas)
1468                  pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+ &
1469                                         zdqcocloud(ig,l,igcm_co_ice)
1470               ENDDO
1471            ENDDO
1472
1473            ! Increment CO ice surface tracer tendency
1474            DO ig=1,ngrid
1475            dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+  &
1476                                     zdqscocloud(ig,igcm_co_ice)
1477            ENDDO
1478
1479            ! update temperature tendancy
1480            DO ig=1,ngrid
1481               DO l=1,nlayer
1482               pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
1483               ENDDO
1484            ENDDO
1485         ELSE
1486         rice_co(:,:)=0 ! initialization needed for callsedim
1487         END IF  ! of IF (carbox)
1488
1489         ! ----------------------------------------
1490         !   VI.1. Microphysics / Aerosol particles
1491         ! ----------------------------------------
1492         ! Call of microphysics
1493         IF (callmufi) THEN
1494
1495            ! Production for microphysics
1496            IF (call_haze_prod_pCH4) THEN
1497               zdqphot_prec(:,:)   = 0.
1498               zdqphot_ch4(:,:)    = 0.
1499               pdqmufi_prod(:,:,:) = 0.
1500               call hazecloud(ngrid,nlayer,nq,ptimestep,                             &
1501                              pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,pdqmufi_prod, &
1502                              zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1503               pdq(:,:,:) = pdq(:,:,:) + pdqmufi_prod(:,:,:) ! Should be updated
1504            ENDIF ! end call_haze_prod_pCH4
1505
1506            pdqmufi(:,:,:) = 0.
1507
1508            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,pdqmufi_prod,pdqmufi)
1509
1510            pdq(:,:,:) = pdq(:,:,:) + pdqmufi(:,:,:)
1511
1512         ELSE
1513            IF (haze) THEN
1514               zdqphot_prec(:,:) = 0.
1515               zdqphot_ch4(:,:)  = 0.
1516               zdqhaze(:,:,:)    = 0.
1517               ! Forcing to a fixed haze profile if haze_proffix
1518               if (haze_proffix.and.i_haze.gt.0.) then
1519                  call haze_prof(ngrid,nlayer,zzlay,pplay,pt,  &
1520                                 reffrad,profmmr)
1521                  zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))/ptimestep
1522               else
1523                  call hazecloud(ngrid,nlayer,nq,ptimestep,            &
1524                     pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, &
1525                     zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
1526               endif
1527               pdq(:,:,:) = pdq(:,:,:) + zdqhaze(:,:,:) ! Should be updated
1528            ENDIF ! end haze
1529
1530            IF (fast.and.fasthaze) THEN
1531               call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_star, &
1532                        mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot,   &
1533                        fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)
1534               DO ig=1,ngrid
1535                  pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+  &
1536                                                zdqprodhaze(ig,igcm_ch4_gas)
1537                  pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+ &
1538                                             zdqprodhaze(ig,igcm_prec_haze)
1539                  pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+ &
1540                                             zdqprodhaze(ig,igcm_haze))
1541                  qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+ &
1542                                                zdqsprodhaze(ig)*ptimestep
1543               ENDDO
1544            ENDIF ! end fast.and.fasthaze
1545
1546            ! Sedimentation.
1547            if (sedimentation) then
1548               zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
1549               zdqssed(1:ngrid,1:nq)  = 0.0
1550               if (oldplutosedim)then
1551                  call callsedim_pluto(ngrid,nlayer,ptimestep,    &
1552                           pplev,zzlev,pt,pdt,rice_ch4,rice_co, &
1553                           pq,pdq,zdqsed,zdqssed,nq,pphi)
1554               else
1555                  call callsedim(ngrid,nlayer,ptimestep,       &
1556                           pplev,zzlev,pt,pdt,pq,pdq,        &
1557                           zdqsed,zdqssed,nq,pphi)
1558               endif
1559               ! Whether it falls as rain or snow depends only on the surface temperature
1560               pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
1561               dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
1562            end if ! end of 'sedimentation'
1563
1564         ENDIF ! end callmufi
1565
1566  ! ---------------
1567  !   VI.2. Updates
1568  ! ---------------
1569
1570         ! Updating Atmospheric Mass and Tracers budgets.
1571         if(mass_redistrib) then
1572
1573            zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
1574            !    (   zdqevap(1:ngrid,1:nlayer)                          &
1575            !    !   + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas)             &
1576            !    !   + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas)             &
1577            !      + dqvaplscale(1:ngrid,1:nlayer) )
1578
1579            do ig = 1, ngrid
1580               zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
1581            enddo
1582
1583            ! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1584            ! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1585            call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1586
1587            call mass_redistribution(ngrid,nlayer,nq,ptimestep,                     &
1588                                     capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
1589                                     pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
1590                                     zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1591
1592            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
1593            dqsurf(1:ngrid,1:nq)       = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
1594            pdt(1:ngrid,1:nlayer)      = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
1595            pdu(1:ngrid,1:nlayer)      = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
1596            pdv(1:ngrid,1:nlayer)      = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
1597            pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
1598            zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
1599
1600         endif
1601
1602!         call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
1603
1604  ! -----------------------------
1605  !   VI.3. Surface Tracer Update
1606  ! -----------------------------
1607
1608         qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
1609
1610      endif! end of if 'tracer'
1611
1612      if (conservn2) then
1613         write(*,*) 'conservn2 tracer'
1614         call testconservmass(ngrid,nlayer,pplev(:,1)+   &
1615            pdpsrf(:)*ptimestep,qsurf(:,1))
1616      endif
1617
1618      DO ig=1,ngrid
1619        flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)- &
1620                                   qsurf1(ig,igcm_n2))/ptimestep
1621        flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
1622        if (methane) then
1623          flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)- &
1624                                   qsurf1(ig,igcm_ch4_ice))/ptimestep
1625          flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
1626        endif
1627        if (carbox) then
1628          flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)-   &
1629                                   qsurf1(ig,igcm_co_ice))/ptimestep
1630          !flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
1631        endif
1632      ENDDO
1633
1634      !! Special source of haze particle !
1635      ! todo: should be placed in haze and use tendency of n2 instead of flusurf
1636      IF (source_haze) THEN
1637         write(*,*) "Source haze not supported yet."
1638         stop
1639            !  call hazesource(ngrid,nlayer,nq,ptimestep,  &
1640            !                 pplev,flusurf,mu0,zdq_source)
1641
1642             DO iq=1, nq
1643               DO l=1,nlayer
1644                 DO ig=1,ngrid
1645                    pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
1646                 ENDDO
1647               ENDDO
1648             ENDDO
1649      ENDIF
1650
1651!------------------------------------------------
1652!   VII. Surface and sub-surface soil temperature
1653!------------------------------------------------
1654
1655!   For diagnostic
1656!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1657      if (.not.fast) then
1658       DO ig=1,ngrid
1659         rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
1660         sensiblehf1(ig)=rho(ig,1)*cpp*(0.4/log(zzlay(ig,1)/z0))**2* &
1661                        (pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1))**0.5*  &
1662                        (tsurf(ig)-pt(ig,1))
1663         if (calldifv) then
1664            sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
1665         end if
1666       ENDDO
1667      endif
1668
1669
1670! VII.1 Increment surface temperature
1671!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1672      tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
1673
1674      ! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
1675      ! Lellouch et al., 2000,2011
1676      IF (tsurfmax) THEN
1677        DO ig=1,ngrid
1678         if (albedo_equivalent(ig).gt.albmin_ch4.and. &
1679                           qsurf(ig,igcm_n2).eq.0.) then
1680              tsurf(ig)=min(tsurf(ig),54.)
1681         endif
1682        ENDDO
1683      ENDIF
1684
1685! VII.2 Compute soil temperatures and subsurface heat flux.
1686!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1687      if (callsoil) then
1688         call soil(ngrid,nsoilmx,.false.,lastcall,therm_inertia,   &
1689                   ptimestep,tsurf,tsoil,capcal,fluxgrd)
1690      endif
1691
1692      ! ! For output :
1693    !   tidat_out(:,:)=0.
1694    !   DO l=1,nsoilmx
1695    !      tidat_out(:,l)=therm_inertia(:,l)
1696    !   ENDDO
1697
1698      ! Test energy conservation
1699      if(enertest)then
1700         call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)
1701         if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
1702      endif
1703
1704
1705
1706!   VII.3 multiply tendencies of cond/subli for paleo loop only in the
1707!       last Pluto year of the simulation
1708!       Year day must be adapted in the startfi for each object
1709!       Paleo uses year_day to calculate the annual mean tendancies
1710!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1711      IF (paleo) then
1712         if (zday.gt.day_ini+ptime0+nday-year_day) then
1713            DO iq=1,nq
1714             DO ig=1,ngrid
1715               qsurfyear(ig,iq)=qsurfyear(ig,iq)+ &
1716                             (qsurf(ig,iq)-qsurf1(ig,iq))  !kg m-2 !ptimestep
1717             ENDDO
1718            ENDDO
1719         endif
1720      endif
1721
1722!   VII.4 Glacial flow at each timestep glastep or at lastcall
1723!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1724      IF (fast.and.glaflow) THEN
1725         if ((mod(zday-day_ini-ptime0,glastep)).lt.1. &
1726                                                   .or.lastcall) then
1727           IF (lastcall) then
1728            dstep=mod(zday-day_ini-ptime0,glastep)*daysec
1729           else
1730            dstep=glastep*daysec
1731           endif
1732           zdqflow(:,:)=qsurf(:,:)
1733           IF (paleo) then
1734             call spreadglacier_paleo(ngrid,nq,qsurf, &
1735                                    phisfinew,dstep,tsurf)
1736           else
1737             call spreadglacier_simple(ngrid,nq,qsurf,dstep)
1738           endif
1739           zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
1740
1741           if (conservn2) then
1742            write(*,*) 'conservn2 glaflow'
1743            call testconservmass(ngrid,nlayer,pplev(:,1)+ &
1744            pdpsrf(:)*ptimestep,qsurf(:,1))
1745           endif
1746
1747         endif
1748      ENDIF
1749
1750!---------------------------------------------------
1751!   VIII. Perform diagnostics and write output files
1752!---------------------------------------------------
1753
1754      ! Note : For output only: the actual model integration is performed in the dynamics.
1755
1756
1757      ! Temperature, zonal and meridional winds.
1758      zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
1759      zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
1760      zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
1761
1762      !! Recast thermal plume vertical velocity array for outputs
1763      !! AF24: removed
1764
1765      ! Diagnostic.
1766      zdtdyn(1:ngrid,1:nlayer)     = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
1767      ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
1768
1769      zdudyn(1:ngrid,1:nlayer)     = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
1770      zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
1771
1772      if(firstcall)then
1773         zdtdyn(1:ngrid,1:nlayer)=0.0
1774         zdudyn(1:ngrid,1:nlayer)=0.0
1775      endif
1776
1777      ! Dynamical heating diagnostic
1778      fluxdyn(:)=0.0
1779      if (.not.fast) then
1780        do ig=1,ngrid
1781           fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp
1782        enddo
1783      endif
1784
1785      ! Tracers.
1786      zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
1787
1788      ! Surface pressure.
1789      ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
1790      call planetwide_sumval(ps(:)*cell_area(:)/totarea_planet,globave)
1791
1792      ! pressure density !pluto specific
1793      IF (.not.fast) then !
1794         do ig=1,ngrid
1795            do l=1,nlayer
1796               zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
1797               zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
1798               rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
1799            enddo
1800            zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
1801         enddo
1802      ENDIF
1803
1804
1805      ! Surface and soil temperature information
1806      call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1807      call planetwide_minval(tsurf(:),Ts2)
1808      call planetwide_maxval(tsurf(:),Ts3)
1809      if(callsoil)then
1810         TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
1811         if (is_master) then
1812            print*,'          ave[Tsurf]             min[Tsurf]             max[Tsurf]             ave[Tdeep]'
1813            print*,Ts1,Ts2,Ts3,TsS
1814         end if
1815      else
1816         if (is_master) then
1817            print*,'          ave[Tsurf]             min[Tsurf]             max[Tsurf]'
1818            print*,Ts1,Ts2,Ts3
1819         endif
1820      end if
1821
1822
1823      ! Check the energy balance of the simulation during the run
1824      if(corrk)then
1825
1826         call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1827         call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1828         call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1829         call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1830         call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1831         do ig=1,ngrid
1832            if(fluxtop_dn(ig).lt.0.0)then
1833               print*,'fluxtop_dn has gone crazy'
1834               print*,'fluxtop_dn=',fluxtop_dn(ig)
1835               print*,'tau_col=',tau_col(ig)
1836               print*,'dtau_aer=',dtau_aer(ig,:,:)
1837               print*,'temp=   ',pt(ig,:)
1838               print*,'pplay=  ',pplay(ig,:)
1839               call abort
1840            endif
1841         end do
1842
1843         if(ngrid.eq.1)then
1844            DYN=0.0
1845         endif
1846
1847         if (is_master) then
1848            print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
1849            print*, ISR,ASR,OLR,GND,DYN
1850         endif
1851
1852         if(enertest .and. is_master)then
1853            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1854            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1855            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1856         endif
1857
1858         if(meanOLR .and. is_master)then
1859            if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1860               ! to record global radiative balance
1861               open(92,file="rad_bal.out",form='formatted',position='append')
1862               write(92,*) zday,ISR,ASR,OLR
1863               close(92)
1864               open(93,file="tem_bal.out",form='formatted',position='append')
1865               if(callsoil)then
1866                         write(93,*) zday,Ts1,Ts2,Ts3,TsS
1867               else
1868                  write(93,*) zday,Ts1,Ts2,Ts3
1869               endif
1870               close(93)
1871            endif
1872         endif
1873
1874      endif ! end of 'corrk'
1875
1876
1877      ! Diagnostic to test radiative-convective timescales in code.
1878      if(testradtimes)then
1879         open(38,file="tau_phys.out",form='formatted',position='append')
1880         ig=1
1881         do l=1,nlayer
1882            write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1883         enddo
1884         close(38)
1885         print*,'As testradtimes enabled,'
1886         print*,'exiting physics on first call'
1887         call abort
1888      endif
1889
1890
1891      ! Compute column amounts (kg m-2) if tracers are enabled.
1892      if(tracer)then
1893         qcol(1:ngrid,1:nq)=0.0
1894         do iq=1,nq
1895            do ig=1,ngrid
1896               qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
1897            enddo
1898         enddo
1899
1900      endif ! end of 'tracer'
1901
1902      if (methane) then
1903         IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
1904           DO ig=1,ngrid
1905             vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)* &
1906                              mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
1907           ENDDO
1908         ELSE
1909           DO ig=1,ngrid
1910 !           compute vmr methane
1911             vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)* &
1912                   g/ps(ig)*mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
1913 !           compute density methane
1914             DO l=1,nlayer
1915                zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
1916             ENDDO
1917           ENDDO
1918         ENDIF
1919      endif
1920
1921      if (carbox) then
1922         IF (fast) then
1923           DO ig=1,ngrid
1924             vmr_co(ig)=zq(ig,1,igcm_co_gas)*   &
1925                              mmol(igcm_n2)/mmol(igcm_co_gas)*100.
1926           ENDDO
1927         ELSE
1928          DO ig=1,ngrid
1929 !           compute vmr CO
1930             vmr_co(ig)=qcol(ig,igcm_co_gas)*   &
1931                   g/ps(ig)*mmol(igcm_n2)/mmol(igcm_co_gas)*100.
1932 !           compute density CO
1933             DO l=1,nlayer
1934                zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
1935             ENDDO
1936          ENDDO
1937         ENDIF
1938      endif
1939
1940      zrho_haze(:,:)=0.
1941      zdqrho_photprec(:,:)=0.
1942      IF (haze.and.optichaze) then
1943         DO ig=1,ngrid
1944            DO l=1,nlayer
1945               zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
1946               zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
1947            ENDDO
1948         ENDDO
1949       ENDIF
1950
1951      IF (fasthaze) then
1952         DO ig=1,ngrid
1953            qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
1954            qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
1955         ENDDO
1956      ENDIF
1957
1958 !     Info about Ls, declin...
1959      IF (fast) THEN
1960         if (is_master) write(*,*),'Ls=',zls*180./pi,' dec=',declin*180./pi
1961         if (is_master) write(*,*),'zday=',zday,' ps=',globave
1962         IF (lastcall) then
1963            if (is_master) write(*,*),'lastcall'
1964         ENDIF
1965      ELSE
1966         if (is_master) write(*,*),'Ls=',zls*180./pi,'decli=',declin*180./pi,'zday=',zday
1967      ENDIF
1968
1969      lecttsoil=0 ! default value for lecttsoil
1970      call getin_p("lecttsoil",lecttsoil)
1971      IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
1972      ! save tsoil temperature profile for 1D profile
1973         OPEN(13,file='proftsoil.out',form='formatted')
1974         DO i=1,nsoilmx
1975            write(13,*) tsoil(1,i)
1976         ENDDO
1977         CLOSE(13)
1978      ENDIF
1979
1980      if (is_master) print*,'--> Ls =',zls*180./pi
1981
1982      if(lastcall) then
1983         IF (grid_type==unstructured) THEN !IF DYNAMICO
1984            ! DYNAMICO: no need to add a dynamics time step to ztime_fin
1985            ztime_fin = ptime
1986         ELSE ! LMDZ
1987            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1988         ENDIF ! of IF (grid_type==unstructured)
1989         !! Update surface ice distribution to iterate to steady state if requested
1990         !! AF24: removed
1991
1992         ! endif
1993         if (paleo) then
1994            ! time range for tendencies of ice flux qsurfyear
1995            zdt_tot=year_day   ! Last year of simulation
1996
1997            masslost(:)=0.
1998            massacc(:)=0.
1999
2000            DO ig=1,ngrid
2001               ! update new reservoir of ice on the surface
2002               DO iq=1,nq
2003                ! kg/m2 to be sublimed or condensed during paleoyears
2004                qsurfyear(ig,iq)=qsurfyear(ig,iq)* &
2005                           paleoyears*365.25/(zdt_tot*daysec/86400.)
2006
2007               ! special case if we sublime the entire reservoir
2008               !! AF: TODO : fix following lines (real_area), using line below:
2009            ! call planetwide_sumval((-qsurfyear(:,iq)-qsurf(:,iq))*cell_area(:),masslost)
2010
2011               !  IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
2012               !    masslost(iq)=masslost(iq)+real_area(ig)*   &
2013               !          (-qsurfyear(ig,iq)-qsurf(ig,iq))
2014               !    qsurfyear(ig,iq)=-qsurf(ig,iq)
2015               !  ENDIF
2016
2017               !  IF (qsurfyear(ig,iq).gt.0.) THEN
2018               !    massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq)
2019               !  ENDIF
2020
2021
2022               ENDDO
2023            ENDDO
2024
2025            DO ig=1,ngrid
2026               DO iq=1,nq
2027                 qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
2028                 IF (qsurfyear(ig,iq).gt.0.) THEN
2029                  qsurfpal(ig,iq)=qsurfpal(ig,iq)- &
2030                       qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
2031                 ENDIF
2032               ENDDO
2033            ENDDO
2034            ! Finally ensure conservation of qsurf
2035            DO iq=1,nq
2036               call planetwide_sumval(qsurf(:,iq)*cell_area(:)/totarea_planet,globaveice(iq))
2037               call planetwide_sumval(qsurfpal(:,iq)*cell_area(:)/totarea_planet,globavenewice(iq))
2038               IF (globavenewice(iq).gt.0.) THEN
2039                  qsurfpal(:,iq)=qsurfpal(:,iq)* &
2040                                   globaveice(iq)/globavenewice(iq)
2041               ENDIF
2042            ENDDO
2043
2044            ! update new geopotential depending on the ice reservoir
2045            phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
2046            !phisfipal(ig)=phisfi(ig)
2047
2048            if (kbo.or.triton) then !  case of Triton : we do not change the orbital parameters
2049               pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
2050               eccpal=1.-periastr/((periastr+apoastr)/2.)    !no change of ecc
2051               peri_daypal=peri_day ! no change
2052               oblipal=obliquit     ! no change
2053               tpalnew=tpal
2054               adjustnew=adjust
2055
2056            else  ! Pluto
2057               ! update new pday and tpal (Myr) to be set in startfi controle
2058               pdaypal=int(day_ini+paleoyears*365.25/6.3872)
2059               tpalnew=tpal+paleoyears*1.e-6  ! Myrs
2060
2061               ! update new N2 ice adjustment (not tested yet on Pluto)
2062               adjustnew=adjust
2063
2064               ! update milankovitch parameters : obliquity,Lsp,ecc
2065               call calcmilank(tpalnew,oblipal,peri_daypal,eccpal)
2066               !peri_daypal=peri_day
2067               !eccpal=0.009
2068            endif
2069
2070            if (is_master) write(*,*) "Paleo peri=",peri_daypal,"  tpal=",tpalnew
2071            if (is_master) write(*,*) "Paleo eccpal=",eccpal,"  tpal=",tpalnew
2072
2073!----------------------------------------------------------------------
2074!        Writing NetCDF file  "RESTARTFI" at the end of the run
2075!----------------------------------------------------------------------
2076!        Note: 'restartfi' is stored just before dynamics are stored
2077!              in 'restart'. Between now and the writing of 'restart',
2078!              there will have been the itau=itau+1 instruction and
2079!              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2080!              thus we store for time=time+dtvr
2081
2082            ! create restartfi
2083            if (ngrid.ne.1) then
2084               print*, "physdem1pal not yet implemented"
2085               stop
2086               !TODO: import this routine from pluto.old
2087               ! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, &
2088               !      ptimestep,pdaypal, &
2089               !      ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, &
2090               !      cell_area,albedodat,therm_inertia,zmea,zstd,zsig, &
2091               !      zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal,  &
2092               !      peri_daypal)
2093            endif
2094         else ! 'paleo'
2095
2096            if (ngrid.ne.1) then
2097               write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
2098
2099               call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
2100                           ptimestep,ztime_fin,tsurf,                &
2101                           tsoil,therm_inertia,emis,albedo,q2,qsurf)
2102            endif
2103
2104         endif ! end of 'paleo'
2105      endif ! end of 'lastcall'
2106
2107!------------------------------------------------------------------------------
2108!           OUTPUT in netcdf file "DIAGFI.NC",
2109!           containing any variable for diagnostic
2110!
2111!             Note 1 : output with  period "ecritphy", set in "run.def"
2112!             Note 2 : writediagfi can also be called from any other subroutine
2113!                      for any variable, but its preferable to keep all the
2114!                      calls in one place ...
2115!------------------------------------------------------------------------------
2116
2117      !-------- General 1D variables
2118
2119      call write_output("Ls","solar longitude","deg",zls*180./pi)
2120      ! call write_output("Lss","sub solar longitude","deg",zlss*180./pi)
2121      call write_output("RA","right ascension","deg",right_ascen*180./pi)
2122      call write_output("Declin","solar declination","deg",declin*180./pi)
2123      call write_output("dist_star","dist_star","AU",dist_star)
2124      call write_output("globave","surf press","Pa",globave)
2125
2126      !-------- General 2D variables
2127
2128      call write_output("tsurf","Surface temperature","K",tsurf)
2129      call write_output("ps","Surface pressure","Pa",ps)
2130      call write_output("emis","Emissivity","",emis)
2131      !if (grid_type == regular_lonlat) then
2132      !     call write_output("area","Mesh area","m2", &
2133      !                      cell_area_for_lonlat_outputs)
2134      !   else ! unstructured grid (e.g. dynamico)
2135      !     call write_output("area","Mesh area","m2",cell_area)
2136      !endif
2137
2138      if (fast) then
2139         call write_output("fluxrad","fluxrad","W m-2",fluxrad)
2140         call write_output("fluxgrd","fluxgrd","W m-2",fluxgrd)
2141         ! call write_output("dplanck","dplanck","W.s m-2 K-1",dplanck)
2142         ! "soil" variables
2143         call write_output("capcal","capcal","W.s m-2 K-1",capcal)
2144         call write_output("tsoil","tsoil","K",tsoil)
2145         call write_output("therm_inertia","therm_inertia","S.I.",therm_inertia)
2146      endif
2147
2148      ! Total energy balance diagnostics
2149      if(callrad)then
2150         call write_output("ALB","Surface albedo"," ",albedo_equivalent)
2151         call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw)
2152         call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn)
2153         call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw)
2154         call write_output("GND","heat flux from ground","W m-2",fluxgrd)
2155         if (.not.fast) then
2156           call write_output("DYN","dynamical heat input","W m-2",fluxdyn)
2157         endif
2158      endif ! end of 'callrad'
2159
2160      !-------- General 3D variables
2161
2162      if (.not.fast) then
2163         if (check_physics_outputs) then
2164            ! Check the validity of updated fields at the end of the physics step
2165            call check_physics_fields("physiq:", zt, zu, zv, pplev, zq)
2166         endif
2167
2168         call write_output("zzlay","Midlayer altitude", "m",zzlay(:,:))
2169         call write_output("zzlev","Interlayer altitude", "m",zzlev(:,1:nlayer))
2170         !call write_output('pphi','Geopotential',' ',pphi)
2171
2172         call write_output("temperature","temperature","K",zt)
2173         call write_output("teta","potential temperature","K",zh)
2174         call write_output("u","Zonal wind","m.s-1",zu)
2175         call write_output("v","Meridional wind","m.s-1",zv)
2176         call write_output("w","Vertical wind","m.s-1",pw)
2177         call write_output("p","Pressure","Pa",pplay)
2178         call write_output("omega","omega","Pa/s",omega)
2179      endif
2180
2181      if(enertest) then
2182         if (calldifv) then
2183            call write_output("q2","turbulent kinetic energy","J.kg^-1",q2)
2184            call write_output("sensibFlux","sensible heat flux","w.m^-2",sensibFlux)
2185         endif
2186         if (corrk) then
2187            call write_output("dEzradsw","radiative heating","w.m^-2",dEzRadsw)
2188            call write_output("dEzradlw","radiative heating","w.m^-2",dEzRadlw)
2189         endif
2190      endif ! end of 'enertest'
2191
2192      ! Diagnostics of optical thickness
2193      ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
2194      if (diagdtau.and..not.oldplutocorrk) then
2195            do nw=1,L_NSPECTV
2196               write(str2,'(i2.2)') nw
2197               call write_output('dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',int_dtauv(:,nlayer:1:-1,nw))
2198            enddo
2199            do nw=1,L_NSPECTI
2200               write(str2,'(i2.2)') nw
2201               call write_output('dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',int_dtaui(:,nlayer:1:-1,nw))
2202            enddo
2203      endif
2204
2205      ! Temporary inclusions for heating diagnostics.
2206      if (.not.fast) then
2207        call write_output("zdtsw","SW heating","K s-1",zdtsw)
2208        call write_output("zdtlw","LW heating","K s-1",zdtlw)
2209        call write_output("dtrad","radiative heating","K s-1",dtrad)
2210        call write_output("zdtdyn","Dyn. heating","K s-1",zdtdyn)
2211        call write_output("zdudyn","Dyn. U","m s-2",zdudyn)
2212        call write_output("zdtconduc","tendancy conduc","K s-1",zdtconduc)
2213        call write_output("zdumolvis","tendancy molvis","m s-1",zdumolvis)
2214        call write_output("zdvmolvis","tendancy molvis","m s-1",zdvmolvis)
2215        call write_output("zdtdif","tendancy T diff","K s-1",zdtdif)
2216        call write_output("zdtsdif","tendancy Ts diff","K s-1",zdtsdif)
2217        call write_output("zdtadj","tendancy T adj","K s-1",zdtadj)
2218      endif
2219
2220      ! Output tracers.
2221      if (tracer) then
2222
2223         do iq=1,nq
2224            if (.not.fast) then
2225              call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq))
2226            endif
2227            call write_output(trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
2228                           'kg m^-2',qcol(:,iq) )
2229            call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
2230                         'kg m^-2',qsurf(:,iq) )
2231         enddo ! end of 'nq' loop
2232
2233         ! N2 cycle
2234         call write_output('n2_iceflux','n2_iceflux',"kg m^-2 s^-1",flusurf(:,igcm_n2) )
2235         if (.not.fast) then
2236            call write_output("zdtc","tendancy T cond N2","K s-1",zdtc)
2237            call write_output("zdtsurfc","tendancy Ts cond N2","K s-1",zdtsurfc)
2238            call write_output("zduc","tendancy U cond N2","m s-1",zduc)
2239            call write_output("zdvc","tendancy V cond N2","m s-1",zdvc)
2240            call write_output("zdqc_n2","tendancy tracer cond N2","kg kg-1 s-1",zdqc(:,:,1))
2241            call write_output("zdqsc_n2","tendancy tracer surf cond N2","kg kg-1 s-1",zdqsc(:,1))
2242            call write_output("zdqdif_n2","tendancy tracer diff","kg kg-1 s-1",zdqdif(:,:,1))
2243            call write_output("zdqsdif_n2","tendancy tracer surf diff","kg kg-1 s-1",zdqsdif(:,1))
2244            call write_output("zdqadj_n2","tendancy tracer adj","K s-1",zdqadj(:,:,1))
2245         endif
2246
2247         ! CH4 cycle
2248         if (methane) then
2249
2250            call write_output('ch4_iceflux','ch4_iceflux',&
2251                              "kg m^-2 s^-1",flusurf(:,igcm_ch4_ice) )
2252            call write_output("vmr_ch4","vmr_ch4","%",vmr_ch4)
2253
2254            if (.not.fast) then
2255               call write_output("zrho_ch4","zrho_ch4","kg.m-3",zrho_ch4(:,:))
2256               !call write_output("rice_ch4","ch4 ice mass mean radius","m",rice_ch4)
2257               !call write_output("zq1temp_ch4"," "," ",zq1temp_ch4)
2258               !call write_output("qsat_ch4"," "," ",qsat_ch4)
2259               !call write_output("qsat_ch4_l1"," "," ",qsat_ch4_l1)
2260
2261               ! 3D Tendancies
2262               call write_output("zdqcn2_ch4","zdq condn2 ch4","",&
2263                           zdqc(:,:,igcm_ch4_gas))
2264               call write_output("zdqdif_ch4","zdqdif ch4","",&
2265                           zdqdif(:,:,igcm_ch4_gas))
2266               call write_output("zdqsdif_ch4_ice","zdqsdif ch4","",&
2267                           zdqsdif(:,igcm_ch4_ice))
2268               call write_output("zdqadj_ch4","zdqadj ch4","",&
2269                           zdqadj(:,:,igcm_ch4_gas))
2270            endif
2271
2272            if (sedimentation) then
2273               call write_output("zdqsed_ch4","zdqsed ch4","",&
2274                              zdqsed(:,:,igcm_ch4_gas))
2275               call write_output("zdqssed_ch4","zdqssed ch4","",&
2276                              zdqssed(:,igcm_ch4_gas))
2277            endif
2278
2279            if (metcloud.and.(.not.fast)) then
2280               call write_output("zdtch4cloud","ch4 cloud","K s-1",&
2281                           zdtch4cloud)
2282               call write_output("zdqch4cloud_gas","ch4 cloud","kg kg-1 s-1",&
2283                           zdqch4cloud(:,:,igcm_ch4_gas))
2284               call write_output("zdqch4cloud_ice","ch4 cloud","kg kg-1 s-1",&
2285                           zdqch4cloud(:,:,igcm_ch4_ice))
2286            endif
2287
2288         endif
2289
2290         ! CO cycle
2291         if (carbox) then
2292            ! call write_output("zdtcocloud","tendancy T cocloud","K",zdtcocloud)
2293            call write_output('co_iceflux','co_iceflux',&
2294                               "kg m^-2 s^-1",flusurf(:,igcm_co_ice) )
2295            call write_output("vmr_co","vmr_co","%",vmr_co)
2296            if (.not.fast) THEN
2297               call write_output("zrho_co","zrho_co","kg.m-3",zrho_co(:,:))
2298            endif
2299         endif
2300
2301         ! Haze
2302         if (haze) then
2303
2304            if (haze_radproffix)then
2305               call write_output('haze_reff','haze_reff','m',reffrad(:,:,1))
2306            end if
2307            !call write_output("zrho_haze","zrho_haze","kg.m-3",zrho_haze(:,:))
2308            !call write_output("zdqhaze_col","zdqhaze col","kg/m2/s",&
2309            !                        zdqhaze_col(:))
2310
2311            ! 3D Tendencies
2312            call write_output("zdqrho_photprec","zdqrho_photprec",&
2313                        "kg.m-3.s-1",zdqrho_photprec(:,:))
2314            call write_output("zdqphot_prec","zdqphot_prec","",&
2315                                                zdqphot_prec(:,:))
2316            call write_output("zdqhaze_ch4","zdqhaze_ch4","",&
2317                     zdqhaze(:,:,igcm_ch4_gas))
2318            call write_output("zdqhaze_prec","zdqhaze_prec","",&
2319                     zdqhaze(:,:,igcm_prec_haze))
2320            call write_output("zdqphot_ch4","zdqphot_ch4","",&
2321                                                zdqphot_ch4(:,:))
2322            call write_output("zdqconv_prec","zdqconv_prec","",&
2323                                                zdqconv_prec(:,:))
2324
2325            if (igcm_haze.ne.0) then
2326               call write_output("zdqhaze_haze","zdqhaze_haze","",&
2327                        zdqhaze(:,:,igcm_haze))
2328               if (sedimentation) then
2329                  call write_output("zdqssed_haze","zdqssed haze",&
2330                     "kg/m2/s",zdqssed(:,igcm_haze))
2331               endif
2332            endif
2333
2334           if (optichaze) then
2335              call write_output("tau_col",&
2336                 "Total aerosol optical depth","opacity",tau_col)
2337           endif
2338
2339         endif ! end haze
2340
2341         if (callmufi) then
2342            ! Tracers:
2343            call write_output("m0as","Spherical mode 0th order moment","m-3",zq(:,:,micro_indx(1))*int2ext(:,:))
2344            call write_output("m3as","Spherical mode 3rd order moment","m3.m-3",zq(:,:,micro_indx(2))*int2ext(:,:))
2345            call write_output("m0af","Fractal mode 0th order moment","m-3",zq(:,:,micro_indx(3))*int2ext(:,:))
2346            call write_output("m3af","Fractal mode 3rd order moment","m3.m-3",zq(:,:,micro_indx(4))*int2ext(:,:))
2347
2348            ! Diagnostics:
2349            call write_output("rcs","Spherical mode characteristic radius","m",mp2m_rc_sph(:,:))
2350            call write_output("rcf","Fractal mode characteristic radius","m",mp2m_rc_fra(:,:))
2351
2352           if (optichaze) then
2353              call write_output("tau_col",&
2354                 "Total aerosol optical depth","opacity",tau_col)
2355           endif ! end optichaze
2356         endif ! end callmufi
2357
2358      endif ! end of 'tracer'
2359
2360      ! Output spectrum.
2361      if(specOLR.and.corrk)then
2362         call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
2363         call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
2364         call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu)
2365      endif
2366
2367! XIOS outputs
2368#ifdef CPP_XIOS
2369      ! Send fields to XIOS: (NB these fields must also be defined as
2370      ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2371      CALL send_xios_field("controle",tab_cntrl_mod,1)
2372
2373      CALL send_xios_field("ap",ap,1)
2374      CALL send_xios_field("bp",bp,1)
2375      CALL send_xios_field("aps",aps,1)
2376      CALL send_xios_field("bps",bps,1)
2377
2378      if (lastcall.and.is_omp_master) then
2379        write(*,*) "physiq: call xios_context_finalize"
2380        call xios_context_finalize
2381      endif
2382#endif
2383
2384      if (check_physics_outputs) then
2385         ! Check the validity of updated fields at the end of the physics step
2386         call check_physics_fields("end of physiq:", zt, zu, zv, pplev, zq)
2387      endif
2388
2389      icount=icount+1
2390
2391      end subroutine physiq
2392
2393   end module physiq_mod
Note: See TracBrowser for help on using the repository browser.