source: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F @ 3726

Last change on this file since 3726 was 3726, checked in by emillour, 3 months ago

Mars PCM:
Turn "callkeys.h" into module "callkeys_mod.F90"
EM

File size: 163.1 KB
Line 
1      MODULE physiq_mod
2
3      IMPLICIT NONE
4
5      CONTAINS
6
7      SUBROUTINE physiq(
8     $            ngrid,nlayer,nq
9     $            ,firstcall,lastcall
10     $            ,pday,ptime,ptimestep
11     $            ,pplev,pplay,pphi
12     $            ,pu,pv,pt,pq
13     $            ,flxw
14     $            ,pdu,pdv,pdt,pdq,pdpsrf)
15
16      use watercloud_mod, only: watercloud, zdqcloud, zdqscloud
17      use calchim_mod, only: calchim, ichemistry, zdqchim, zdqschim
18      use watersat_mod, only: watersat
19      use co2condens_mod, only: co2condens, CO2cond_ps
20      use co2cloud_mod, only: co2cloud
21      use callradite_mod, only: callradite
22      use callsedim_mod, only: callsedim
23      use rocketduststorm_mod, only: rocketduststorm, dustliftday
24      use calcstormfract_mod, only: calcstormfract
25      use topmons_mod, only: topmons,topmons_setup
26      use nltecool_mod, only: nltecool
27      use nlte_tcool_mod, only: nlte_tcool
28      use blendrad_mod, only: blendrad
29      use nlthermeq_mod, only: nlthermeq
30      use thermosphere_mod, only: thermosphere
31      use param_read_e107_mod, only: param_read_e107
32      use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
33     &                      igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
34     &                      igcm_hdo_vap, igcm_hdo_ice,
35     &                      igcm_ccn_mass, igcm_ccn_number,
36     &                      igcm_ccnco2_mass, igcm_ccnco2_number,
37     &                      igcm_ccnco2_h2o_mass_ice,
38     &                      igcm_ccnco2_h2o_mass_ccn,
39     &                      igcm_ccnco2_h2o_number,
40     &                      igcm_ccnco2_meteor_mass,
41     &                      igcm_ccnco2_meteor_number,
42     &                      igcm_dust_mass, igcm_dust_number, igcm_h2o2,
43     &                      nuice_ref, rho_ice, rho_dust, ref_r0,
44     &                      igcm_he, igcm_stormdust_mass,
45     &                      igcm_stormdust_number, igcm_topdust_mass,
46     &                      igcm_topdust_number,
47     &                      qperemin
48      use comsoil_h, only: inertiedat, inertiesoil,! dat: soil thermal inertia for present climate, inertiesoil is the TI read in the start
49     &                     tsoil, nsoilmx,!number of subsurface layers
50     &                     mlayer,layer, ! soil mid layer depths
51     &                     nqsoil,qsoil  ! adsorption
52      use geometry_mod, only: longitude, latitude, cell_area,
53     &                        cell_area_for_lonlat_outputs,longitude_deg
54      use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
55      use surfdat_h, only: phisfi, albedodat, z0, albedo_h2o_cap,
56     &                     albedo_h2o_frost, frost_albedo_threshold,
57     &                     frost_metam_threshold, tsurf, emis, capcal,
58     &                     fluxgrd, qsurf, watercap, watercaptag,
59     &                     perennial_co2ice
60      use comsaison_h, only: dist_sol, declin, zls,
61     &                       mu0, fract, local_time
62      use solarlong_mod, only: solarlong
63      use nirdata_mod, only: NIR_leedat
64      use nirco2abs_mod, only: nirco2abs
65      use surfacearea_mod, only: surfacearea
66      use slope_mod, only: theta_sl, psi_sl, getslopes, param_slope
67      use conc_mod, only: init_r_cp_mu, update_r_cp_mu_ak, rnew,
68     &                    cpnew, mmean
69      use time_phylmdz_mod, only: steps_per_sol
70      use time_phylmdz_mod, only: iphysiq, ecritstart, daysec
71      use dimradmars_mod, only: aerosol, totcloudfrac,
72     &                          dtrad, fluxrad_sky, fluxrad, albedo,
73     &                          naerkind, iaer_dust_doubleq,
74     &                          iaer_stormdust_doubleq, iaer_h2o_ice,
75     &                          flux_1AU
76      use dust_param_mod, only: doubleq, lifting, callddevil,
77     &                          tauscaling, odpref, dustbin,
78     &                          dustscaling_mode, dust_rad_adjust,
79     &                          freedust, reff_driven_IRtoVIS_scenario
80      use dustdevil_mod, only: dustdevil
81      use turb_mod, only: q2, wstar, ustar, sensibFlux,
82     &                    zmax_th, hfmax_th, turb_resolved
83      use planete_h, only: aphelie, periheli, year_day, peri_day,
84     &                     obliquit
85      use planete_h, only: iniorbit
86      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
87      USE calldrag_noro_mod, ONLY: calldrag_noro
88      USE vdifc_mod, ONLY: vdifc
89      use param_v4_h, only: nreact,n_avog,
90     &                      fill_data_thermos, allocate_param_thermos
91      use iono_h, only: allocate_param_iono
92      use compute_dtau_mod, only: compute_dtau
93      use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
94      use nonoro_gwd_mix_mod, only: nonoro_gwd_mix, calljliu_gwimix
95      use check_fields_mod, only: check_physics_fields
96      use surfini_mod, only: surfini
97#ifdef MESOSCALE
98      use comsoil_h, only: mlayer,layer
99      use surfdat_h, only: z0_default
100      use comm_wrf
101#else
102      USE planetwide_mod, ONLY: planetwide_maxval, planetwide_minval,
103     &                          planetwide_sumval
104      use phyredem, only: physdem0, physdem1
105      use phyetat0_mod, only: phyetat0, tab_cntrl_mod
106      use wstats_mod, only: callstats, wstats, mkstats
107      use eofdump_mod, only: eofdump
108      USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
109      USE mod_phys_lmdz_omp_data, ONLY: is_omp_master
110      USE time_phylmdz_mod, ONLY: day_end
111#endif
112
113#ifdef CPP_XIOS
114      use xios_output_mod, only: initialize_xios_output,
115     &                           update_xios_timestep,
116     &                           send_xios_field
117      use wxios, only: wxios_context_init, xios_context_finalize
118#endif
119      USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured,
120     &                             regular_lonlat
121      use ioipsl_getin_p_mod, only: getin_p
122      use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
123     &                        subslope_dist,iflat,sky_slope,
124     &                        major_slope,compute_meshgridavg,
125     &                        ini_comslope_h
126      use write_output_mod, only: write_output
127      use pbl_parameters_mod, only: pbl_parameters
128      use lmdz_atke_turbulence_ini, only : atke_ini
129      use waterice_tifeedback_mod, only : waterice_tifeedback
130      use callkeys_mod, only: calladj, calltherm, callatke, calldifv
131      use callkeys_mod, only: callrichsl, tke_heat_flux
132      use callkeys_mod, only: calllott, calllott_nonoro, calleofdump
133      use callkeys_mod, only: callrad, callnlte, callnirco2, nircorr
134      use callkeys_mod, only: diurnal, season, iradia, nltemodel
135      use callkeys_mod, only: water, activice, microphys, CLFvarying
136      use callkeys_mod, only: hdo, co2clouds, co2useh2o, meteo_flux
137      use callkeys_mod, only: callsoil, callslope, callcond
138      use callkeys_mod, only: tituscap, surfaceice_tifeedback
139      use callkeys_mod, only: refill_watercap, poreice_tifeedback
140      use callkeys_mod, only: cst_cap_albedo
141      use callkeys_mod, only: rdstorm, dustinjection
142      use callkeys_mod, only: topflows, dustiropacity
143      use callkeys_mod, only: sedimentation, scavenging
144      use callkeys_mod, only: photochem, callthermos
145      use callkeys_mod, only: startphy_file
146
147      IMPLICIT NONE
148c=======================================================================
149c
150c   subject:
151c   --------
152c
153c   Organisation of the physical parametrisations of the LMD
154c   martian atmospheric general circulation model.
155c
156c   The GCM can be run without or with tracer transport
157c   depending on the value of Logical "tracer" in file  "callphys.def"
158c   Tracers may be water vapor, ice OR chemical species OR dust particles
159c
160c   SEE comments in initracer.F about numbering of tracer species...
161c
162c   It includes:
163c
164c      1. Initialization:
165c      1.1 First call initializations
166c      1.2 Initialization for every call to physiq
167c      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
168c      2. Compute radiative transfer tendencies
169c         (longwave and shortwave) for CO2 and aerosols.
170c      3. Gravity wave and subgrid scale topography drag :
171c      4. Vertical diffusion (turbulent mixing):
172c      5. Convective adjustment
173c      6. Condensation and sublimation of carbon dioxide.
174c      7.  TRACERS :
175c       7a. water, water ice, co2 ice (clouds)
176c       7b. call for photochemistry when tracers are chemical species
177c       7c. other scheme for tracer (dust) transport (lifting, sedimentation)
178c       7d. updates (CO2 pressure variations, surface budget)
179c      8. Contribution to tendencies due to thermosphere
180c      9. Surface and sub-surface temperature calculations
181c     10. Write outputs :
182c           - "startfi", "histfi" (if it's time)
183c           - Saving statistics (if "callstats = .true.")
184c           - Dumping eof (if "calleofdump = .true.")
185c           - Output any needed variables in "diagfi"
186c     11. Diagnostic: mass conservation of tracers
187c
188c   author:
189c   -------
190c           Frederic Hourdin    15/10/93
191c           Francois Forget     1994
192c           Christophe Hourdin  02/1997
193c           Subroutine completly rewritten by F.Forget (01/2000)
194c           Introduction of the photochemical module: S. Lebonnois (11/2002)
195c           Introduction of the thermosphere module: M. Angelats i Coll (2002)
196c           Water ice clouds: Franck Montmessin (update 06/2003)
197c           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
198c             Nb: See callradite.F for more information.
199c           Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
200c           jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
201c
202c           10/16 J. Audouard: modifications for CO2 clouds scheme
203
204c   arguments:
205c   ----------
206c
207c   input:
208c   ------
209c    ecri                  period (in dynamical timestep) to write output
210c    ngrid                 Size of the horizontal grid.
211c                          All internal loops are performed on that grid.
212c    nlayer                Number of vertical layers.
213c    nq                    Number of advected fields
214c    firstcall             True at the first call
215c    lastcall              True at the last call
216c    pday                  Number of days counted from the North. Spring
217c                          equinoxe.
218c    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
219c    ptimestep             timestep (s)
220c    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
221c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
222c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
223c    pu(ngrid,nlayer)      u component of the wind (ms-1)
224c    pv(ngrid,nlayer)      v component of the wind (ms-1)
225c    pt(ngrid,nlayer)      Temperature (K)
226c    pq(ngrid,nlayer,nq)   Advected fields
227c    pudyn(ngrid,nlayer)    |
228c    pvdyn(ngrid,nlayer)    | Dynamical temporal derivative for the
229c    ptdyn(ngrid,nlayer)    | corresponding variables
230c    pqdyn(ngrid,nlayer,nq) |
231c    flxw(ngrid,nlayer)      vertical mass flux (kg/s) at layer lower boundary
232c
233c   output:
234c   -------
235c
236c    pdu(ngrid,nlayer)       |
237c    pdv(ngrid,nlayer)       |  Temporal derivative of the corresponding
238c    pdt(ngrid,nlayer)       |  variables due to physical processes.
239c    pdq(ngrid,nlayer,nq)    |
240c    pdpsrf(ngrid)           |
241
242c
243c=======================================================================
244c
245c    0.  Declarations :
246c    ------------------
247
248      include "netcdf.inc"
249
250c Arguments :
251c -----------
252
253c   inputs:
254c   -------
255      INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns
256      INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers
257      INTEGER,INTENT(in) :: nq ! number of tracers
258      LOGICAL,INTENT(in) :: firstcall ! signals first call to physics
259      LOGICAL,INTENT(in) :: lastcall ! signals last call to physics
260      REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0
261      REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
262      REAL,INTENT(in) :: ptimestep ! physics timestep (s)
263      REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
264      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
265      REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
266      REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
267      REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
268      REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
269      REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
270      REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s) at lower boundary of layer
271
272c   outputs:
273c   --------
274c     physical tendencies
275      REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
276      REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
277      REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
278      REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
279      REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
280
281c Local saved variables:
282c ----------------------
283      INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
284      INTEGER,SAVE :: icount  ! Counter of calls to physiq during the run.
285      REAL,SAVE :: time_phys
286
287!$OMP THREADPRIVATE(day_ini,icount,time_phys)
288
289#ifdef DUSTSTORM
290      REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
291#endif
292
293c     Variables used by the water ice microphysical scheme:
294      REAL rice(ngrid,nlayer)    ! Water ice geometric mean radius (m)
295      REAL nuice(ngrid,nlayer)   ! Estimated effective variance
296                                     !   of the size distribution
297      real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
298      real rhocloud(ngrid,nlayer)  ! Cloud density (kg.m-3)
299      real rsedcloudco2(ngrid,nlayer) ! CO2 Cloud sedimentation radius
300      real rhocloudco2(ngrid,nlayer)  ! CO2 Cloud density (kg.m-3)
301      real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the
302                                  !   size distribution
303      REAL inertiesoil_tifeedback(ngrid,nsoilmx,nslope) ! Time varying subsurface
304                                      ! thermal inertia (J.s-1/2.m-2.K-1)
305                                      ! (used only when tifeedback surface or pore =.true.)
306c     Variables used by the CO2 clouds microphysical scheme:
307      DOUBLE PRECISION riceco2(ngrid,nlayer)   ! co2 ice geometric mean radius (m)
308      real zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
309      real zdqssed_ccn(ngrid,nq)  ! CCN flux at the surface  (kg.m-2.s-1)
310      real, dimension(ngrid,nlayer) :: zcondicea_co2microp
311c     Variables used by the photochemistry
312      REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
313      REAL surfice(ngrid,nlayer)  !  ice surface area (m2/m3, if photochemistry)
314c     Variables used by the slope model
315      REAL sl_lct, sl_lat
316      REAL sl_tau, sl_alb, sl_the, sl_psi
317      REAL sl_fl0, sl_flu
318      REAL sl_ra, sl_di0
319      REAL sky
320      REAL fluxsurf_dir_dn_sw(ngrid) ! Incident direct solar flux on Mars at surface (W.m-2)
321
322      REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant
323
324c Local variables :
325c -----------------
326      REAL CBRT
327      EXTERNAL CBRT
328
329!      CHARACTER*80 fichier
330      INTEGER l,ig,ierr,igout,iq,isoil
331
332      REAL fluxsurf_lw(ngrid,nslope)      !incident LW (IR) surface flux (W.m-2)
333      REAL fluxsurf_dn_sw(ngrid,2,nslope) ! Incident SW (solar) surface flux (W.m-2)
334      REAL fluxsurf_up_sw(ngrid,2) ! Reflected SW (solar) surface flux (W.m-2)
335      REAL fluxtop_lw(ngrid)       !Outgoing LW (IR) flux to space (W.m-2)
336      REAL fluxtop_dn_sw(ngrid,2) ! Incoming SW (solar) flux from space (W.m-2)
337      REAL fluxtop_up_sw(ngrid,2) ! Outgoing SW (solar) flux to space (W.m-2)
338      REAL tau_pref_scenario(ngrid) ! prescribed dust column visible opacity
339                                    ! at odpref
340      REAL IRtoVIScoef(ngrid) ! conversion coefficient to apply on
341                              ! scenario absorption IR (9.3um) CDOD
342                              ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
343      REAL tau_pref_gcm(ngrid) ! dust column visible opacity at odpref in the GCM
344c     rocket dust storm
345      REAL totstormfract(ngrid)     ! fraction of the mesh where the dust storm is contained
346      logical clearatm              ! clearatm used to calculate twice the radiative
347                                    ! transfer when rdstorm is active :
348                                    !            - in a mesh with stormdust and background dust (false)
349                                    !            - in a mesh with background dust only (true)
350c     entrainment by mountain top dust flows
351      logical nohmons               ! nohmons used to calculate twice the radiative
352                                    ! transfer when topflows is active :
353                                    !            - in a mesh with topdust and background dust (false)
354                                    !            - in a mesh with background dust only (true)
355
356      REAL tau(ngrid,naerkind)     ! Column dust optical depth at each point
357                                   ! AS: TBD: this one should be in a module !
358      REAL zday                      ! date (time since Ls=0, in martian days)
359      REAL zzlay(ngrid,nlayer)     ! altitude at the middle of the layers
360      REAL zzlev(ngrid,nlayer+1)   ! altitude at layer boundaries
361      REAL gz(ngrid,nlayer)        ! variation of g with altitude from aeroid surface
362!      REAL latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
363
364c     Tendancies due to various processes:
365      REAL dqsurf(ngrid,nq,nslope)  ! tendency for tracers on surface (Kg/m2/s)
366      REAL zdtlw(ngrid,nlayer)     ! (K/s)
367      REAL zdtsw(ngrid,nlayer)     ! (K/s)
368      REAL pdqrds(ngrid,nlayer,nq) ! tendency for dust after rocketduststorm
369
370      REAL zdtnirco2(ngrid,nlayer) ! (K/s)
371      REAL zdtnlte(ngrid,nlayer)   ! (K/s)
372      REAL zdtsurf(ngrid,nslope)            ! (K/s)
373      REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer)
374      REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer)  ! (m.s-2)
375      REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid,nslope)         ! (K/s)
376      REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer)  ! (m.s-2)
377      REAL zdhadj(ngrid,nlayer)                           ! (K/s)
378      REAL zdtgw(ngrid,nlayer)                            ! (K/s)
379      REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer)    ! (m.s-2)
380      REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid,nslope)
381      REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
382
383      REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq,nslope)
384      REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
385      REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
386      REAL zdqadj(ngrid,nlayer,nq)
387      REAL zdqc(ngrid,nlayer,nq)
388      REAL zdqcloudco2(ngrid,nlayer,nq)
389      REAL zdqsc(ngrid,nq,nslope)
390
391      REAL zdteuv(ngrid,nlayer)    ! (K/s)
392      REAL zdtconduc(ngrid,nlayer) ! (K/s)
393      REAL zdumolvis(ngrid,nlayer)
394      REAL zdvmolvis(ngrid,nlayer)
395      real zdqmoldiff(ngrid,nlayer,nq)
396      real*8 PhiEscH,PhiEscH2,PhiEscD
397
398      REAL dwatercap(ngrid,nslope), dwatercap_dif(ngrid,nslope)     ! (kg/m-2)
399
400c     Local variable for local intermediate calcul:
401      REAL zflubid(ngrid,nslope)
402      REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
403      REAL zdum1(ngrid,nlayer)
404      REAL zdum2(ngrid,nlayer)
405      REAL ztim1,ztim2,ztim3, z1,z2
406      REAL ztime_fin
407      REAL zdh(ngrid,nlayer)
408      REAL zh(ngrid,nlayer)      ! potential temperature (K)
409      REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
410      INTEGER length
411      PARAMETER (length=100)
412      REAL tlaymean ! temporary value of mean layer temperature for zzlay
413
414c     Variables for the total dust for diagnostics
415      REAL qdusttotal(ngrid,nlayer) !it equals to dust + stormdust
416
417c local variables only used for diagnostic (output in file "diagfi" or "stats")
418c -----------------------------------------------------------------------------
419      REAL ps(ngrid), zt(ngrid,nlayer)
420      REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
421      REAL zq(ngrid,nlayer,nq)
422
423      REAL fluxtop_dn_sw_tot(ngrid), fluxtop_up_sw_tot(ngrid)
424      REAL fluxsurf_dn_sw_tot(ngrid,nslope), fluxsurf_up_sw_tot(ngrid)
425      character*2 str2
426!      character*5 str5
427      real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
428      real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
429      real rstormdust(ngrid,nlayer) ! stormdust geometric mean radius (m)
430      real rtopdust(ngrid,nlayer)   ! topdust geometric mean radius (m)
431      integer igmin, lmin
432
433      ! pplev and pplay are dynamical inputs and must not be modified in the physics.
434      ! instead, use zplay and zplev :
435      REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
436!      REAL zstress(ngrid),cd
437      real rho(ngrid,nlayer)  ! density
438      real vmr(ngrid,nlayer)  ! volume mixing ratio
439      real rhopart(ngrid,nlayer) ! number density of a given species
440      real colden(ngrid,nq)     ! vertical column of tracers
441      real mass(nq)             ! global mass of tracers (g)
442      REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
443      REAL mstormdtot(ngrid)    ! Total mass of stormdust tracer (kg/m2)
444      REAL mdusttot(ngrid)      ! Total mass of dust tracer (kg/m2)
445      REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
446      REAL mtotco2(ngrid)      ! Total mass of co2, including ice at the surface (kg/m2)
447      REAL vaptotco2(ngrid)     ! Total mass of co2 vapor (kg/m2)
448      REAL icetotco2(ngrid)     ! Total mass of co2 ice (kg/m2)
449      REAL Nccntot(ngrid)       ! Total number of ccn (nbr/m2)
450      REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
451      REAL rave(ngrid)          ! Mean water ice effective radius (m)
452      REAL opTES(ngrid,nlayer)  ! abs optical depth at 825 cm-1
453      REAL tauTES(ngrid)        ! column optical depth at 825 cm-1
454      REAL Qabsice                ! Water ice absorption coefficient
455      REAL taucloudtes(ngrid)  ! Cloud opacity at infrared
456                               !   reference wavelength using
457                               !   Qabs instead of Qext
458                               !   (direct comparison with TES)
459      REAL mtotD(ngrid)          ! Total mass of HDO vapor (kg/m2)
460      REAL icetotD(ngrid)        ! Total mass of HDO ice (kg/m2)
461      REAL DoH_vap(ngrid,nlayer) !D/H ratio
462      REAL DoH_ice(ngrid,nlayer) !D/H ratio
463      REAL DoH_surf(ngrid) !D/H ratio surface
464
465      REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
466      REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
467      REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
468      REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
469      REAL nccn(ngrid,nlayer)  ! true n ccn (kg/kg)
470      REAL qccn(ngrid,nlayer)  ! true q ccn (kg/kg)
471c     definition tendancies of stormdust tracers
472      REAL rdsdqdustsurf(ngrid) ! surface q stormdust flux (kg/m2/s)
473      REAL rdsdndustsurf(ngrid) ! surface n stormdust flux (number/m2/s)
474      REAL rdsndust(ngrid,nlayer) ! true n stormdust (kg/kg)
475      REAL rdsqdust(ngrid,nlayer) ! true q stormdust (kg/kg)
476      REAL wspeed(ngrid,nlayer+1) ! vertical velocity stormdust tracer
477      REAL wtop(ngrid,nlayer+1) ! vertical velocity topdust tracer
478      REAL dsodust(ngrid,nlayer) ! density scaled opacity for background dust
479      REAL dsords(ngrid,nlayer) ! density scaled opacity for stormdust
480      REAL dsotop(ngrid,nlayer) ! density scaled opacity for topdust
481
482c Test 1d/3d scavenging
483      REAL satu(ngrid,nlayer)  ! satu ratio for output
484      REAL zqsat(ngrid,nlayer) ! saturation
485
486! Added for new NLTE scheme
487      real co2vmr_gcm(ngrid,nlayer)
488      real n2vmr_gcm(ngrid,nlayer)
489      real ovmr_gcm(ngrid,nlayer)
490      real covmr_gcm(ngrid,nlayer)
491      integer ierr_nlte
492      real*8  varerr
493
494c  Non-oro GW tendencies
495      REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
496      REAL d_t_hin(ngrid,nlayer)
497      REAL d_u_mix(ngrid,nlayer), d_v_mix(ngrid,nlayer)
498      REAL d_t_mix(ngrid,nlayer), zdq_mix(ngrid,nlayer,nq)
499
500c  Diagnostics 2D of gw_nonoro
501      REAL zustrhi(ngrid), zvstrhi(ngrid)
502c Variables for PBL
503      REAL zz1(ngrid)
504      REAL lmax_th_out(ngrid)
505      REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
506      REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
507      INTEGER lmax_th(ngrid),n_out,n
508      CHARACTER(50) zstring
509      REAL dtke_th(ngrid,nlayer+1)
510      REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
511      REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
512      REAL u_out1(ngrid)
513      REAL T_out1(ngrid)
514      REAL, ALLOCATABLE, DIMENSION(:) :: z_out     ! height of interpolation between z0 and z1 [meters]
515      REAL tstar(ngrid)  ! friction velocity and friction potential temp
516      REAL vhf(ngrid), vvv(ngrid)
517      real qdustrds0(ngrid,nlayer),qdustrds1(ngrid,nlayer)
518      real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer)
519      real qdusttotal0(ngrid),qdusttotal1(ngrid)
520
521c sub-grid scale water ice clouds (A. Pottier 2013)
522      logical clearsky
523      ! flux for the part without clouds
524      real zdtswclf(ngrid,nlayer)
525      real zdtlwclf(ngrid,nlayer)
526      real fluxsurf_lwclf(ngrid)
527      real fluxsurf_dn_swclf(ngrid,2),fluxsurf_up_swclf(ngrid,2)
528      real fluxtop_lwclf(ngrid)
529      real fluxtop_dn_swclf(ngrid,2),fluxtop_up_swclf(ngrid,2)
530      real taucloudtesclf(ngrid)
531      real tf_clf, ntf_clf ! tf: fraction of clouds, ntf: fraction without clouds
532      real rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m)
533C test de conservation de la masse de CO2
534      REAL co2totA
535      REAL co2totB
536      REAL co2conservation
537
538c entrainment by mountain top dust flows above sub-grid scale topography
539      REAL pdqtop(ngrid,nlayer,nq) ! tendency for dust after topmons
540
541c when no startfi file is asked for init
542      real alpha,lay1 ! coefficients for building layers
543      integer iloop
544
545      ! flags to trigger extra sanity checks
546      logical,save :: check_physics_inputs=.false.
547      logical,save :: check_physics_outputs=.false.
548
549!$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs)
550
551c Sub-grid scale slopes
552      real :: tsurf_meshavg(ngrid) ! Surface temperature grid box averaged [K]
553      real :: albedo_meshavg(ngrid,2) ! albedo temperature grid box averaged [1]
554      real :: emis_meshavg(ngrid,2) ! emis temperature grid box averaged [1]
555      real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2]
556      real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie
557      integer :: islope
558      real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting
559
560      logical :: write_restart
561
562! Variable for ice table
563      REAL :: rhowater_surf(ngrid,nslope)            ! Water density at the surface [kg/m^3]
564      REAL :: rhowater_surf_sat(ngrid,nslope)        ! Water density at the surface at saturation [kg/m^3]
565      REAL :: rhowater_soil(ngrid,nsoilmx,nslope)    ! Water density in soil layers [kg/m^3]
566      REAL,PARAMETER :: alpha_clap_h2o = 28.9074     ! Coeff for Clapeyron law [/]
567      REAL,PARAMETER :: beta_clap_h2o = -6143.7      ! Coeff for Clapeyron law [K]
568      REAL :: pvap_surf(ngrid)                       ! Water vapor partial pressure in first layer [Pa]
569      REAL,PARAMETER :: m_co2 = 44.01E-3             ! CO2 molecular mass [kg/mol]
570      REAL,PARAMETER :: m_noco2 = 33.37E-3           ! Non condensible mol mass [kg/mol]
571      REAL :: ztmp1,ztmp2                            ! intermediate variables to compute the mean molar mass of the layer
572      REAL :: pore_icefraction(ngrid,nsoilmx,nslope) ! ice filling fraction in the pores
573! Variable for the computation of the TKE with parameterization from ATKE working group
574      REAL :: viscom                              ! kinematic molecular viscosity for momentum
575      REAL :: viscoh                              ! kinematic molecular viscosity for heat
576
577c=======================================================================
578      pdq(:,:,:) = 0.
579
580c 1. Initialisation:
581c -----------------
582c  1.1   Initialisation only at first call
583c  ---------------------------------------
584      IF (firstcall) THEN
585
586         call getin_p("check_physics_inputs",check_physics_inputs)
587         call getin_p("check_physics_outputs",check_physics_outputs)
588
589c        variables set to 0
590c        ~~~~~~~~~~~~~~~~~~
591         aerosol(:,:,:)=0
592         dtrad(:,:)=0
593
594#ifndef MESOSCALE
595         fluxrad(:,:)=0
596         wstar(:)=0.
597#endif
598
599#ifdef CPP_XIOS
600         ! Initialize XIOS context
601         write(*,*) "physiq: call wxios_context_init"
602         CALL wxios_context_init
603#endif
604
605c        read startfi
606c        ~~~~~~~~~~~~
607#ifndef MESOSCALE
608
609! GCM. Read netcdf initial physical parameters.
610         CALL phyetat0 ("startfi.nc",0,0,
611     &         nsoilmx,ngrid,nlayer,nq,nqsoil,
612     &         day_ini,time_phys,
613     &         tsurf,tsoil,albedo,emis,
614     &         q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,
615     &         watercap,perennial_co2ice,
616     &         def_slope,def_slope_mean,subslope_dist)
617
618!     Sky view:
619       DO islope=1,nslope
620       sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
621       END DO
622!     Determine the 'flatest' slopes
623       iflat = 1
624       DO islope=2,nslope
625         IF(abs(def_slope_mean(islope)).lt.
626     &      abs(def_slope_mean(iflat)))THEN
627           iflat = islope
628         ENDIF
629       ENDDO
630       write(*,*)'Flat slope for islope = ',iflat
631       write(*,*)'corresponding criterium = ',def_slope_mean(iflat)
632
633#else
634! MESOSCALE. Supposedly everything is already set in modules.
635! So we just check. And we fill day_ini
636      write(*,*)"check: --- in physiq.F"
637      write(*,*)"check: rad,cpp,g,r,rcp,daysec"
638      write(*,*)rad,cpp,g,r,rcp,daysec
639      write(*,*)'check: tsurf ',tsurf(1,:),tsurf(ngrid,:)
640      write(*,*)'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:)
641      write(*,*)'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
642      write(*,*)'check: midlayer,layer ', mlayer(:),layer(:)
643      write(*,*)'check: tracernames ', noms
644      write(*,*)'check: emis ',emis(1,:),emis(ngrid,:)
645      write(*,*)'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
646      write(*,*)'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:)
647      write(*,*)'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:)
648      !!!
649      day_ini = pday
650      !!! a couple initializations (dummy for mesoscale) done in phyetat0
651      !!! --- maybe this should be done in update_inputs_physiq_mod
652
653      tauscaling(:)=1.0 !! probably important
654      totcloudfrac(:)=1.0
655      DO islope = 1,nslope
656      albedo(:,1,islope)=albedodat(:)
657      albedo(:,2,islope)=albedo(:,1,islope)
658      inertiesoil(:,:,islope) = inertiedat(:,:)
659      watercap(:,:)=0.0
660      ENDDO
661#endif
662#ifndef MESOSCALE
663         if (.not.startphy_file) then
664           ! starting without startfi.nc and with callsoil
665           ! is not yet possible as soildepth default is not defined
666           if (callsoil) then
667              ! default mlayer distribution, following a power law:
668              !  mlayer(k)=lay1*alpha**(k-1/2)
669              lay1=2.e-4
670              alpha=2
671              do iloop=0,nsoilmx-1
672              mlayer(iloop)=lay1*(alpha**(iloop-0.5))
673          enddo
674              lay1=sqrt(mlayer(0)*mlayer(1))
675              alpha=mlayer(1)/mlayer(0)
676              do iloop=1,nsoilmx
677                 layer(iloop)=lay1*(alpha**(iloop-1))
678              enddo
679           endif
680           ! additionnal "academic" initialization of physics
681           do islope = 1,nslope
682             tsurf(:,islope)=pt(:,1)
683           enddo
684           write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!"
685           do isoil=1,nsoilmx
686             tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:)
687           enddo
688           write(*,*) "Physiq: initializing inertiedat !!"
689           inertiedat(:,:)=400.
690           inertiesoil(:,:,:)=400.
691           write(*,*) "Physiq: initializing day_ini to pdat !"
692           day_ini=pday
693        endif
694#endif
695         if (pday.ne.day_ini) then
696           write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
697     &                "physics and dynamics"
698           write(*,*) "dynamics day [pday]: ",pday
699           write(*,*) "physics day [day_ini]: ",day_ini
700           call abort_physic("physiq","dynamics day /= physics day",1)
701         endif
702
703         write (*,*) 'In physiq day_ini =', day_ini
704
705c        initialize tracers
706c        ~~~~~~~~~~~~~~~~~~
707         CALL initracer(ngrid,nq,qsurf)
708
709c        Initialize albedo and orbital calculation
710c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
711         CALL surfini(ngrid,nslope,qsurf)
712         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
713c        initialize soil
714c        ~~~~~~~~~~~~~~~
715         IF (callsoil) THEN
716c           Thermal inertia feedback:
717            IF (surfaceice_tifeedback.or.poreice_tifeedback) THEN
718             DO islope = 1,nslope
719                CALL waterice_tifeedback(ngrid,nsoilmx,nslope,
720     s               qsurf(:,igcm_h2o_ice,:),pore_icefraction,
721     s               inertiesoil_tifeedback)
722             ENDDO
723                CALL soil(ngrid,nsoilmx,firstcall,
724     s             inertiesoil_tifeedback,
725     s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
726            ELSE
727                CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
728     s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
729            ENDIF ! of IF (tifeedback)
730         ELSE
731            write(*,*)
732     &     'PHYSIQ WARNING! Thermal conduction in the soil turned off'
733            DO ig=1,ngrid
734               capcal(ig,:)=1.e5
735               fluxgrd(ig,:)=0.
736            ENDDO
737         ENDIF
738         icount=1
739
740#ifndef MESOSCALE
741c        Initialize thermospheric parameters
742c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
743
744         if (callthermos) then
745            call fill_data_thermos
746            call allocate_param_thermos(nlayer)
747            call allocate_param_iono(nlayer,nreact)
748            call param_read_e107
749         endif
750#endif
751
752c        Initialize rnew cpnew and mmean as constant
753c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
754         call init_r_cp_mu(ngrid,nlayer)
755
756         if(callnlte.and.nltemodel.eq.2) call nlte_setup
757         if(callnirco2.and.nircorr.eq.1) call NIR_leedat
758
759
760        IF (water.AND.(ngrid.NE.1)) THEN
761          write(*,*)"physiq: water_param Surface water frost albedo:",
762     .                  albedo_h2o_frost
763          write(*,*)"physiq: water_param Surface watercap albedo:",
764     .                  albedo_h2o_cap
765        ENDIF
766
767#ifndef MESOSCALE
768      ! no need to compute slopes when in 1D; it is an input
769      if (ngrid /= 1 .and. callslope) call getslopes(ngrid,phisfi)
770      if (ecritstart.GT.0) then
771             call physdem0("restartfi.nc",longitude,latitude,
772     &                   nsoilmx,ngrid,nlayer,nq,
773     &                   ptimestep,pday,0.,cell_area,
774     &                   albedodat,inertiedat,def_slope,
775     &                   subslope_dist)
776      else
777             call physdem0("restartfi.nc",longitude,latitude,
778     &                   nsoilmx,ngrid,nlayer,nq,
779     &                   ptimestep,float(day_end),0.,cell_area,
780     &                   albedodat,inertiedat,def_slope,
781     &                   subslope_dist)
782      endif
783
784c        Initialize mountain mesh fraction for the entrainment by top flows param.
785c        ~~~~~~~~~~~~~~~
786        if (topflows) call topmons_setup(ngrid)
787
788c        Parameterization of the ATKE routine
789c        ~~~~~~~~~~~~~~~
790        if (callatke) then
791           viscom = 0.001
792           viscoh = 0.001
793           CALL atke_ini(g, r, pi, cpp, 0., viscom, viscoh)
794        endif
795
796#endif
797
798#ifdef CPP_XIOS
799        ! XIOS outputs
800        write(*,*) "physiq firstcall: call initialize_xios_output"
801        call initialize_xios_output(pday,ptime,ptimestep,daysec,
802     &                              presnivs,pseudoalt,mlayer)
803#endif
804      ENDIF        !  (end of "if firstcall")
805
806      if (check_physics_inputs) then
807        ! Check the validity of input fields coming from the dynamics
808        call check_physics_fields("begin physiq:",pt,pu,pv,pplev,pq)
809      endif
810
811c ---------------------------------------------------
812c 1.2   Initializations done at every physical timestep:
813c ---------------------------------------------------
814c
815#ifdef CPP_XIOS
816      ! update XIOS time/calendar
817      call update_xios_timestep
818#endif
819
820c     Initialize various variables
821c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
822      pdv(:,:)=0
823      pdu(:,:)=0
824      pdt(:,:)=0
825      pdq(:,:,:)=0
826      pdpsrf(:)=0
827      zflubid(:,:)=0
828      zdtsurf(:,:)=0
829      dqsurf(:,:,:)=0
830      dsodust(:,:)=0.
831      dsords(:,:)=0.
832      dsotop(:,:)=0.
833      dwatercap(:,:)=0
834
835      call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
836     &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
837
838      ! Dust scenario conversion coefficient from IRabs to VISext
839      IRtoVIScoef(1:ngrid)=2.6 ! initialized with former value from Montabone et al 2015
840                               ! recomputed in aeropacity if reff_driven_IRtoVIS_scenario=.true.
841
842#ifdef DUSTSTORM
843      pq_tmp(:,:,:)=0
844#endif
845      igout=ngrid/2+1
846
847
848      zday=pday+ptime ! compute time, in sols (and fraction thereof)
849      ! Compute local time at each grid point
850      DO ig=1,ngrid
851       local_time(ig)=modulo(1.+(zday-INT(zday))
852     &                +(longitude_deg(ig)/15)/24,1.)
853      ENDDO
854c     Compute Solar Longitude (Ls) :
855c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
856      if (season) then
857         call solarlong(zday,zls)
858      else
859         call solarlong(float(day_ini),zls)
860      end if
861
862c     Initialize pressure levels
863c     ~~~~~~~~~~~~~~~~~~
864      zplev(:,:) = pplev(:,:)
865      zplay(:,:) = pplay(:,:)
866      ps(:) = pplev(:,1)
867
868#ifndef MESOSCALE
869c-----------------------------------------------------------------------
870c    1.2.1 Compute mean mass, cp, and R
871c    update_r_cp_mu_ak outputs rnew(ngrid,nlayer), cpnew(ngrid,nlayer)
872c   , mmean(ngrid,nlayer) and Akknew(ngrid,nlayer)
873c    --------------------------------
874
875      if(photochem.or.callthermos) then
876         call update_r_cp_mu_ak(ngrid,nlayer,nq,
877     &                       zplay,pt,pdt,pq,pdq,ptimestep)
878      endif
879#endif
880
881c     Compute geopotential at interlayers
882c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
883c     ponderation des altitudes au niveau des couches en dp/p
884cc ------------------------------------------
885    !Calculation zzlev & zzlay for first layer
886      DO ig=1,ngrid
887         zzlay(ig,1)=-(log(pplay(ig,1)/ps(ig)))*rnew(ig,1)*pt(ig,1)/g
888         zzlev(ig,1)=0
889         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
890         gz(ig,1)=g
891
892        DO l=2,nlayer
893         ! compute "mean" temperature of the layer
894          if(pt(ig,l) .eq. pt(ig,l-1)) then
895             tlaymean=pt(ig,l)
896          else
897             tlaymean=(pt(ig,l)- pt(ig,l-1))/log(pt(ig,l)/pt(ig,l-1))
898          endif
899
900          ! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
901          gz(ig,l)=g*(rad**2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))**2
902
903          zzlay(ig,l)=zzlay(ig,l-1)-
904     &   (log(pplay(ig,l)/pplay(ig,l-1))*rnew(ig,l)*tlaymean/gz(ig,l))
905
906          z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
907          z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
908          zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
909        ENDDO
910      ENDDO
911
912!     Potential temperature calculation not the same in physiq and dynamic
913
914c     Compute potential temperature
915c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
916      DO l=1,nlayer
917         DO ig=1,ngrid
918            zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
919            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
920         ENDDO
921      ENDDO
922
923
924      ! Compute vertical velocity (m/s) from vertical mass flux
925      ! w = F / (rho*area) and rho = P/(r*T)
926      ! but first linearly interpolate mass flux to mid-layers
927      do l=1,nlayer-1
928       pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
929      enddo
930      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
931      do l=1,nlayer
932       pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /
933     &               (pplay(1:ngrid,l)*cell_area(1:ngrid))
934       ! NB: here we use r and not rnew since this diagnostic comes
935       ! from the dynamics
936      enddo
937
938      ! test for co2 conservation with co2 microphysics
939      if (igcm_co2_ice.ne.0) then
940        ! calculates the amount of co2 at the beginning of physics
941        co2totA = 0.
942        do ig=1,ngrid
943          do l=1,nlayer
944             co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
945     &             (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
946     &        +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
947          end do
948          do islope = 1,nslope
949          co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
950     &    subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
951          enddo
952        end do
953      else
954        co2totA = 0.
955        do ig=1,ngrid
956          do l=1,nlayer
957             co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
958     &             (pq(ig,l,igcm_co2)
959     &        +pdq(ig,l,igcm_co2)*ptimestep)
960          end do
961          do islope = 1,nslope
962          co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
963     &    subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
964          enddo
965        end do
966      endif ! of if (igcm_co2_ice.ne.0)
967c-----------------------------------------------------------------------
968c    2. Compute radiative tendencies :
969c------------------------------------
970
971      IF (callrad) THEN
972
973c       Local Solar zenith angle
974c       ~~~~~~~~~~~~~~~~~~~~~~~~
975
976        CALL orbite(zls,dist_sol,declin)
977
978        IF (diurnal) THEN
979            ztim1=SIN(declin)
980            ztim2=COS(declin)*COS(2.*pi*(zday-.5))
981            ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
982
983            CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
984     &                  ztim1,ztim2,ztim3, mu0,fract)
985
986        ELSE
987            CALL mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad)
988        ENDIF ! of IF (diurnal)
989
990         IF( MOD(icount-1,iradia).EQ.0) THEN
991
992c          NLTE cooling from CO2 emission
993c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
994           IF(callnlte) then
995              if(nltemodel.eq.0.or.nltemodel.eq.1) then
996                CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
997              else if(nltemodel.eq.2) then
998                co2vmr_gcm(1:ngrid,1:nlayer)=
999     &                      pq(1:ngrid,1:nlayer,igcm_co2)*
1000     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
1001                n2vmr_gcm(1:ngrid,1:nlayer)=
1002     &                      pq(1:ngrid,1:nlayer,igcm_n2)*
1003     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
1004                covmr_gcm(1:ngrid,1:nlayer)=
1005     &                      pq(1:ngrid,1:nlayer,igcm_co)*
1006     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
1007                ovmr_gcm(1:ngrid,1:nlayer)=
1008     &                      pq(1:ngrid,1:nlayer,igcm_o)*
1009     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
1010
1011                 CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
1012     $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1013     $                ovmr_gcm,  zdtnlte,ierr_nlte,varerr )
1014                 if(ierr_nlte.gt.0) then
1015                    write(*,*)
1016     $                'WARNING: nlte_tcool output with error message',
1017     $                'ierr_nlte=',ierr_nlte,'varerr=',varerr
1018                    write(*,*)'I will continue anyway'
1019                 endif
1020
1021                 zdtnlte(1:ngrid,1:nlayer)=
1022     &                             zdtnlte(1:ngrid,1:nlayer)/86400.
1023              endif
1024           ELSE
1025              zdtnlte(:,:)=0.
1026           ENDIF !end callnlte
1027
1028!          Find number of layers for LTE radiation calculations
1029!          (done only once per sol)
1030           IF(MOD((icount-1),steps_per_sol).EQ.0)
1031     &          CALL nlthermeq(ngrid,nlayer,zplev,zplay)
1032
1033c          rocketstorm : compute dust storm mesh fraction
1034           IF (rdstorm) THEN
1035              CALL calcstormfract(ngrid,nlayer,nq,pq,
1036     &                 totstormfract)
1037           ENDIF
1038
1039c          Note: Dustopacity.F has been transferred to callradite.F
1040
1041#ifdef DUSTSTORM
1042!! specific case: save the quantity of dust before adding perturbation
1043
1044       if (firstcall) then
1045        pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass)
1046        pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number)
1047       endif
1048#endif
1049
1050c          Call main radiative transfer scheme
1051c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1052c          Transfer through CO2 (except NIR CO2 absorption)
1053c            and aerosols (dust and water ice)
1054           ! callradite for background dust (out of the rdstorm fraction)
1055           clearatm=.true.
1056           !! callradite for background dust (out of the topflows fraction)
1057           nohmons=.true.
1058
1059c          Radiative transfer
1060c          ------------------
1061           ! callradite for the part with clouds
1062           clearsky=.false. ! part with clouds for both cases CLFvarying true/false
1063           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
1064     &     albedo_meshavg,emis_meshavg,
1065     &     mu0,zplev,zplay,pt,tsurf_meshavg,fract,dist_sol,igout,
1066     &     zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat),
1067     &     fluxsurf_up_sw,
1068     &     fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw,
1069     &     tau_pref_scenario,tau_pref_gcm,
1070     &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
1071     &     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,
1072     &     qsurf_meshavg(:,igcm_co2),rstormdust,rtopdust,totstormfract,
1073     &     clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac)
1074
1075           DO islope=1,nslope
1076             fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat)
1077             fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat)
1078           ENDDO
1079
1080           ! case of sub-grid water ice clouds: callradite for the clear case
1081            IF (CLFvarying) THEN
1082               ! ---> PROBLEMS WITH ALLOCATED ARRAYS
1083               ! (temporary solution in callcorrk: do not deallocate
1084               ! if
1085               ! CLFvarying ...) ?? AP ??
1086               clearsky=.true.
1087               CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
1088     &              albedo_meshavg,emis_meshavg,mu0,zplev,zplay,pt,
1089     &              tsurf_meshavg,fract,
1090     &              dist_sol,igout,zdtlwclf,zdtswclf,
1091     &              fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf,
1092     &              fluxtop_lwclf,fluxtop_dn_swclf,fluxtop_up_swclf,
1093     &              tau_pref_scenario,tau_pref_gcm,tau,aerosol,
1094     &              dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
1095     &              taucloudtesclf,rdust,
1096     &              rice,nuice,riceco2, nuiceco2,
1097     &              qsurf_meshavg(:,igcm_co2),
1098     &              rstormdust,rtopdust,totstormfract,
1099     &              clearatm,dsords,dsotop,
1100     &              nohmons,clearsky,totcloudfrac)
1101               clearsky = .false.  ! just in case.
1102               ! Sum the fluxes and heating rates from cloudy/clear
1103               ! cases
1104               DO ig=1,ngrid
1105                  tf_clf=totcloudfrac(ig)
1106                  ntf_clf=1.-tf_clf
1107                  DO islope=1,nslope
1108                  fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig)
1109     &                          + tf_clf*fluxsurf_lw(ig,islope)
1110                  fluxsurf_dn_sw(ig,1:2,islope) =
1111     &                           ntf_clf*fluxsurf_dn_swclf(ig,1:2)
1112     &                          + tf_clf*fluxsurf_dn_sw(ig,1:2,islope)
1113                  ENDDO
1114                  fluxsurf_up_sw(ig,1:2) =
1115     &                           ntf_clf*fluxsurf_up_swclf(ig,1:2)
1116     &                          + tf_clf*fluxsurf_up_sw(ig,1:2)
1117                  fluxtop_lw(ig)  = ntf_clf*fluxtop_lwclf(ig)
1118     &                                      + tf_clf*fluxtop_lw(ig)
1119                  fluxtop_dn_sw(ig,1:2)=ntf_clf*fluxtop_dn_swclf(ig,1:2)
1120     &                                 + tf_clf*fluxtop_dn_sw(ig,1:2)
1121                  fluxtop_up_sw(ig,1:2)=ntf_clf*fluxtop_up_swclf(ig,1:2)
1122     &                                 + tf_clf*fluxtop_up_sw(ig,1:2)
1123                  taucloudtes(ig) = ntf_clf*taucloudtesclf(ig)
1124     &                                      + tf_clf*taucloudtes(ig)
1125                  zdtlw(ig,1:nlayer) = ntf_clf*zdtlwclf(ig,1:nlayer)
1126     &                                      + tf_clf*zdtlw(ig,1:nlayer)
1127                  zdtsw(ig,1:nlayer) = ntf_clf*zdtswclf(ig,1:nlayer)
1128     &                                      + tf_clf*zdtsw(ig,1:nlayer)
1129               ENDDO
1130
1131            ENDIF ! (CLFvarying)
1132
1133!============================================================================
1134
1135#ifdef DUSTSTORM
1136!! specific case: compute the added quantity of dust for perturbation
1137       if (firstcall) then
1138         pdq(1:ngrid,1:nlayer,igcm_dust_mass)=
1139     &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)
1140     &      - pq_tmp(1:ngrid,1:nlayer,1)
1141     &      + pq(1:ngrid,1:nlayer,igcm_dust_mass)
1142         pdq(1:ngrid,1:nlayer,igcm_dust_number)=
1143     &      pdq(1:ngrid,1:nlayer,igcm_dust_number)
1144     &      - pq_tmp(1:ngrid,1:nlayer,2)
1145     &      + pq(1:ngrid,1:nlayer,igcm_dust_number)
1146       endif
1147#endif
1148
1149c          Outputs for basic check (middle of domain)
1150c          ------------------------------------------
1151           write(*,'("Ls =",f11.6," check lat =",f10.6,
1152     &               " lon =",f11.6)')
1153     &           zls*180./pi,latitude(igout)*180/pi,
1154     &                       longitude(igout)*180/pi
1155
1156           write(*,'(" tau_pref_gcm(",f4.0," Pa) =",f9.6,
1157     &             " tau(",f4.0," Pa) =",f9.6)')
1158     &            odpref,tau_pref_gcm(igout),
1159     &            odpref,tau(igout,1)*odpref/zplev(igout,1)
1160
1161
1162c          ---------------------------------------------------------
1163c          Call slope parameterization for direct and scattered flux
1164c          ---------------------------------------------------------
1165           IF(callslope) THEN
1166           ! assume that in this case, nslope = 1
1167            if(nslope.ne.1) then
1168              call abort_physic(
1169     &      "physiq","callslope=true but nslope.ne.1",1)
1170            endif
1171            write(*,*) 'Slope scheme is on and computing...'
1172            DO ig=1,ngrid
1173              sl_the = theta_sl(ig)
1174              IF (sl_the .ne. 0.) THEN
1175                ztim1=fluxsurf_dn_sw(ig,1,iflat)
1176     &               +fluxsurf_dn_sw(ig,2,iflat)
1177                DO l=1,2
1178                 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
1179                 sl_ra  = pi*(1.0-sl_lct/12.)
1180                 sl_lat = 180.*latitude(ig)/pi
1181                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
1182                 sl_alb = albedo(ig,l,iflat)
1183                 sl_psi = psi_sl(ig)
1184                 sl_fl0 = fluxsurf_dn_sw(ig,l,iflat)
1185                 sl_di0 = 0.
1186                 if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then
1187                  sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
1188                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
1189                  sl_di0 = sl_di0/ztim1
1190                  sl_di0 = fluxsurf_dn_sw(ig,l,iflat)*sl_di0
1191                 endif
1192                 ! you never know (roundup concern...)
1193                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
1194                 !!!!!!!!!!!!!!!!!!!!!!!!!!
1195                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
1196     &                             sl_tau, sl_alb, sl_the, sl_psi,
1197     &                             sl_di0, sl_fl0, sl_flu )
1198                 !!!!!!!!!!!!!!!!!!!!!!!!!!
1199                 fluxsurf_dn_sw(ig,l,1) = sl_flu
1200                ENDDO
1201              !!! compute correction on IR flux as well
1202              sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
1203              fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky
1204              ENDIF
1205            ENDDO
1206          ELSE        ! not calling subslope, nslope might be > 1
1207          DO islope = 1,nslope
1208            sl_the=abs(def_slope_mean(islope))
1209            IF (sl_the .gt. 1e-6) THEN
1210              IF(def_slope_mean(islope).ge.0.) THEN
1211                psi_sl(:) = 0. !Northward slope
1212              ELSE
1213                psi_sl(:) = 180. !Southward slope
1214              ENDIF
1215              DO ig=1,ngrid
1216                ztim1=fluxsurf_dn_sw(ig,1,islope)
1217     s                +fluxsurf_dn_sw(ig,2,islope)
1218                DO l=1,2
1219                 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
1220                 sl_ra  = pi*(1.0-sl_lct/12.)
1221                 sl_lat = 180.*latitude(ig)/pi
1222                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
1223                 sl_alb = albedo(ig,l,islope)
1224                 sl_psi = psi_sl(ig)
1225                 sl_fl0 = fluxsurf_dn_sw(ig,l,islope)
1226                 sl_di0 = 0.
1227                 if (mu0(ig) .gt. 0.) then
1228                  sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
1229                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
1230                  sl_di0 = sl_di0/ztim1
1231                  sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0
1232                 endif
1233                 ! you never know (roundup concern...)
1234                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
1235                 !!!!!!!!!!!!!!!!!!!!!!!!!!
1236                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
1237     &                             sl_tau, sl_alb, sl_the, sl_psi,
1238     &                             sl_di0, sl_fl0, sl_flu )
1239                 !!!!!!!!!!!!!!!!!!!!!!!!!!
1240                 fluxsurf_dn_sw(ig,l,islope) = sl_flu
1241                ENDDO
1242              !!! compute correction on IR flux as well
1243
1244              fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope)
1245     &                                *sky_slope(islope)
1246              ENDDO
1247            ENDIF ! sub grid is not flat
1248          ENDDO ! islope = 1,nslope
1249          ENDIF ! callslope
1250
1251c          CO2 near infrared absorption
1252c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1253           zdtnirco2(:,:)=0
1254           if (callnirco2) then
1255              call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
1256     .                       mu0,fract,declin, zdtnirco2)
1257           endif
1258
1259c          Radiative flux from the sky absorbed by the surface (W.m-2)
1260c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1261           DO ig=1,ngrid
1262             DO islope = 1,nslope
1263               fluxrad_sky(ig,islope) =
1264     $           emis(ig,islope)*fluxsurf_lw(ig,islope)
1265     $         +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope))
1266     $         +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope))
1267             ENDDO
1268           ENDDO
1269
1270c          Net atmospheric radiative heating rate (K.s-1)
1271c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1272           IF(callnlte) THEN
1273              CALL blendrad(ngrid, nlayer, zplay,
1274     &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
1275           ELSE
1276              DO l=1,nlayer
1277                 DO ig=1,ngrid
1278                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
1279     &                          +zdtnirco2(ig,l)
1280                  ENDDO
1281              ENDDO
1282           ENDIF
1283
1284        ENDIF ! of if(mod(icount-1,iradia).eq.0)
1285
1286c    Transformation of the radiative tendencies:
1287c    -------------------------------------------
1288
1289c          Net radiative surface flux (W.m-2)
1290c          ~~~~~~~~~~~~~~~~~~~~~~~~~~
1291
1292c
1293           DO ig=1,ngrid
1294             DO islope = 1,nslope
1295               zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope)
1296               zplanck(ig)=emis(ig,islope)*
1297     $         stephan*zplanck(ig)*zplanck(ig)
1298               fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig)
1299               IF(callslope) THEN
1300                 sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
1301                 fluxrad(ig,nslope)=fluxrad(ig,nslope)+
1302     $                             (1.-sky)*zplanck(ig)
1303               ELSE
1304                 fluxrad(ig,islope)=fluxrad(ig,islope) +
1305     $              (1.-sky_slope(iflat))*emis(ig,iflat)*
1306     $              stephan*tsurf(ig,iflat)**4
1307               ENDIF
1308           ENDDO
1309         ENDDO
1310
1311         DO l=1,nlayer
1312            DO ig=1,ngrid
1313               pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
1314            ENDDO
1315         ENDDO
1316
1317      ENDIF ! of IF (callrad)
1318
1319c     3.1 Rocket dust storm
1320c    -------------------------------------------
1321      IF (rdstorm) THEN
1322         clearatm=.false.
1323         pdqrds(:,:,:)=0.
1324         qdusttotal0(:)=0.
1325         qdusttotal1(:)=0.
1326         do ig=1,ngrid
1327           do l=1,nlayer
1328            zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ! updated potential
1329                                             ! temperature tendency
1330            ! for diagnostics
1331!            qdustrds0(ig,l)=pq(ig,l,igcm_dust_mass)+
1332!     &           pdq(ig,l,igcm_dust_mass)*ptimestep
1333!            qstormrds0(ig,l)=pq(ig,l,igcm_stormdust_mass)+
1334!     &           pdq(ig,l,igcm_stormdust_mass)*ptimestep
1335!            qdusttotal0(ig)=qdusttotal0(ig)+(qdustrds0(ig,l)+
1336!     &          qstormrds0(ig,l))*(zplev(ig,l)-
1337!     &           zplev(ig,l+1))/g
1338           enddo
1339         enddo
1340!         call write_output('qdustrds0','qdust before rds',
1341!     &           'kg/kg ',qdustrds0(:,:))
1342!         call write_output('qstormrds0','qstorm before rds',
1343!     &           'kg/kg ',qstormrds0(:,:))
1344
1345         CALL rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep,
1346     &                      pq,pdq,pt,pdt,zplev,zplay,zzlev,
1347     &                      zzlay,zdtsw,zdtlw,
1348c               for radiative transfer
1349     &                      clearatm,icount,zday,zls,
1350     &                      tsurf_meshavg,qsurf_meshavg(:,igcm_co2),
1351     &                      igout,totstormfract,tauscaling,
1352     &                      dust_rad_adjust,IRtoVIScoef,
1353     &                      albedo_meshavg,emis_meshavg,
1354c               input sub-grid scale cloud
1355     &                      clearsky,totcloudfrac,
1356c               input sub-grid scale topography
1357     &                      nohmons,
1358c               output
1359     &                      pdqrds,wspeed,dsodust,dsords,dsotop,
1360     &                      tau_pref_scenario,tau_pref_gcm)
1361
1362c      update the tendencies of both dust after vertical transport
1363         DO l=1,nlayer
1364          DO ig=1,ngrid
1365           pdq(ig,l,igcm_stormdust_mass)=
1366     &              pdq(ig,l,igcm_stormdust_mass)+
1367     &                      pdqrds(ig,l,igcm_stormdust_mass)
1368           pdq(ig,l,igcm_stormdust_number)=
1369     &              pdq(ig,l,igcm_stormdust_number)+
1370     &                      pdqrds(ig,l,igcm_stormdust_number)
1371
1372           pdq(ig,l,igcm_dust_mass)=
1373     &       pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass)
1374           pdq(ig,l,igcm_dust_number)=
1375     &           pdq(ig,l,igcm_dust_number)+
1376     &                  pdqrds(ig,l,igcm_dust_number)
1377
1378          ENDDO
1379         ENDDO
1380         do l=1,nlayer
1381           do ig=1,ngrid
1382             qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+
1383     &           pdq(ig,l,igcm_dust_mass)*ptimestep
1384             qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+
1385     &           pdq(ig,l,igcm_stormdust_mass)*ptimestep
1386             qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+
1387     &          qstormrds1(ig,l))*(zplev(ig,l)-
1388     &           zplev(ig,l+1))/g
1389           enddo
1390         enddo
1391
1392c        for diagnostics
1393!         call write_output('qdustrds1','qdust after rds',
1394!     &           'kg/kg ',qdustrds1(:,:))
1395!         call write_output('qstormrds1','qstorm after rds',
1396!     &           'kg/kg ',qstormrds1(:,:))
1397!
1398!         call write_output('qdusttotal0','q sum before rds',
1399!     &           'kg/m2 ',qdusttotal0(:))
1400!         call write_output('qdusttotal1','q sum after rds',
1401!     &           'kg/m2 ',qdusttotal1(:))
1402
1403      ENDIF ! end of if(rdstorm)
1404
1405c     3.2 Dust entrained from the PBL up to the top of sub-grid scale topography
1406c    -------------------------------------------
1407      IF (topflows) THEN
1408         clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains
1409         nohmons=.false.
1410         pdqtop(:,:,:)=0.
1411         CALL topmons(ngrid,nlayer,nq,ptime,ptimestep,
1412     &                pq,pdq,pt,pdt,zplev,zplay,zzlev,
1413     &                zzlay,zdtsw,zdtlw,
1414     &                icount,zday,zls,tsurf(:,iflat),
1415     &                qsurf_meshavg(:,igcm_co2),
1416     &                igout,aerosol,tauscaling,dust_rad_adjust,
1417     &                IRtoVIScoef,albedo_meshavg,emis_meshavg,
1418     &                totstormfract,clearatm,
1419     &                clearsky,totcloudfrac,
1420     &                nohmons,
1421     &                pdqtop,wtop,dsodust,dsords,dsotop,
1422     &                tau_pref_scenario,tau_pref_gcm)
1423
1424c      update the tendencies of both dust after vertical transport
1425         DO l=1,nlayer
1426          DO ig=1,ngrid
1427           pdq(ig,l,igcm_topdust_mass)=
1428     & pdq(ig,l,igcm_topdust_mass)+
1429     &                      pdqtop(ig,l,igcm_topdust_mass)
1430           pdq(ig,l,igcm_topdust_number)=
1431     & pdq(ig,l,igcm_topdust_number)+
1432     &                      pdqtop(ig,l,igcm_topdust_number)
1433           pdq(ig,l,igcm_dust_mass)=
1434     & pdq(ig,l,igcm_dust_mass)+ pdqtop(ig,l,igcm_dust_mass)
1435           pdq(ig,l,igcm_dust_number)=
1436     & pdq(ig,l,igcm_dust_number)+pdqtop(ig,l,igcm_dust_number)
1437
1438          ENDDO
1439         ENDDO
1440
1441      ENDIF ! end of if (topflows)
1442
1443c     3.3 Dust injection from the surface
1444c    -------------------------------------------
1445           if (dustinjection.gt.0) then
1446
1447             CALL compute_dtau(ngrid,nlayer,
1448     &                  zday,pplev,tau_pref_gcm,
1449     &                  ptimestep,local_time,IRtoVIScoef,
1450     &                          dustliftday)
1451           endif ! end of if (dustinjection.gt.0)
1452
1453c-----------------------------------------------------------------------
1454c    4. Gravity wave and subgrid scale topography drag :
1455c    -------------------------------------------------
1456
1457
1458      IF(calllott)THEN
1459        CALL calldrag_noro(ngrid,nlayer,ptimestep,
1460     &                 zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
1461
1462        DO l=1,nlayer
1463          DO ig=1,ngrid
1464            pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
1465            pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
1466            pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
1467          ENDDO
1468        ENDDO
1469      ENDIF
1470
1471c-----------------------------------------------------------------------
1472c    5. Vertical diffusion (turbulent mixing):
1473c    -----------------------------------------
1474
1475      IF (calldifv) THEN
1476         DO ig=1,ngrid
1477           DO islope = 1,nslope
1478            zflubid(ig,islope)=fluxrad(ig,islope)
1479     &                        +fluxgrd(ig,islope)
1480           ENDDO
1481         ENDDO
1482         zdum1(:,:)=0
1483         zdum2(:,:)=0
1484         do l=1,nlayer
1485            do ig=1,ngrid
1486               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1487            enddo
1488         enddo
1489
1490c ----------------------
1491c Treatment of a special case : using new surface layer (Richardson based)
1492c without using the thermals in gcm and mesoscale can yield problems in
1493c weakly unstable situations when winds are near to 0. For those cases, we add
1494c a unit subgrid gustiness. Remember that thermals should be used we using the
1495c Richardson based surface layer model.
1496        IF ( .not.calltherm
1497     .       .and. callrichsl
1498     .       .and. .not.turb_resolved) THEN
1499
1500          DO ig=1, ngrid
1501             IF (zh(ig,1) .lt. tsurf_meshavg(ig)) THEN
1502               wstar(ig)=1.
1503               hfmax_th(ig)=0.2
1504             ELSE
1505               wstar(ig)=0.
1506               hfmax_th(ig)=0.
1507             ENDIF
1508          ENDDO
1509        ENDIF
1510
1511c ----------------------
1512
1513         IF (tke_heat_flux .ne. 0.) THEN
1514
1515             zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)*
1516     &                 (-alog(zplay(:,1)/zplev(:,1)))
1517             pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
1518         ENDIF
1519
1520c        Calling vdif (Martian version WITH CO2 condensation)
1521         dwatercap_dif(:,:) = 0.
1522         CALL vdifc(ngrid,nlayer,nsoilmx,nq,nqsoil,zpopsk,
1523     $        ptimestep,capcal,
1524     $        zplay,zplev,zzlay,zzlev,z0,
1525     $        pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,
1526     $        qsoil,pore_icefraction,
1527     $        zdum1,zdum2,zdh,pdq,zflubid,
1528     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
1529     &        zdqdif,zdqsdif,wstar,hfmax_th,
1530     &        zcondicea_co2microp,sensibFlux,
1531     &        dustliftday,local_time,watercap,dwatercap_dif)
1532
1533          DO ig=1,ngrid
1534             zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:)
1535             dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:)
1536          ENDDO
1537
1538          call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
1539     &   albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
1540         IF (.not.turb_resolved) THEN
1541          DO l=1,nlayer
1542            DO ig=1,ngrid
1543               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
1544               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
1545               pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
1546
1547               zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
1548            ENDDO
1549          ENDDO
1550
1551           DO iq=1, nq
1552            DO l=1,nlayer
1553              DO ig=1,ngrid
1554                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
1555              ENDDO
1556            ENDDO
1557           ENDDO
1558           DO iq=1, nq
1559              DO ig=1,ngrid
1560                 dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
1561              ENDDO
1562           ENDDO
1563
1564         ELSE
1565           write (*,*) '******************************************'
1566           write (*,*) '** LES mode: the difv part is only used to'
1567           write (*,*) '**  - provide HFX and UST to the dynamics'
1568           write (*,*) '**  - update TSURF'
1569           write (*,*) '******************************************'
1570           !! Specific treatment for lifting in turbulent-resolving mode (AC)
1571           IF (lifting .and. doubleq) THEN
1572             !! lifted dust is injected in the first layer.
1573             !! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
1574             !! => lifted dust is not incremented before the sedimentation step.
1575             zdqdif(1:ngrid,1,1:nq)=0.
1576             zdqdif(1:ngrid,1,igcm_dust_number) =
1577     .        -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_number)
1578             zdqdif(1:ngrid,1,igcm_dust_mass) =
1579     .        -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_mass)
1580             zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
1581             DO iq=1, nq
1582               IF ((iq .ne. igcm_dust_mass)
1583     &          .and. (iq .ne. igcm_dust_number)) THEN
1584                 zdqsdif(:,iq,:)=0.
1585               ENDIF
1586             ENDDO
1587           ELSE
1588             zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
1589             zdqsdif(1:ngrid,1:nq,1:nslope) = 0.
1590           ENDIF
1591         ENDIF
1592      ELSE
1593         DO ig=1,ngrid
1594           DO islope=1,nslope
1595             zdtsurf(ig,islope)=zdtsurf(ig,islope)+
1596     s   (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope)
1597           ENDDO
1598         ENDDO
1599
1600         IF (turb_resolved) THEN
1601            write(*,*) 'Turbulent-resolving mode !'
1602            write(*,*) 'Please set calldifv to T in callphys.def'
1603            call abort_physic("physiq","turbulent-resolving mode",1)
1604         ENDIF
1605      ENDIF ! of IF (calldifv)
1606
1607c-----------------------------------------------------------------------
1608c   6. Thermals :
1609c   -----------------------------
1610
1611      if(calltherm .and. .not.turb_resolved) then
1612
1613        call calltherm_interface(ngrid,nlayer,nq,igcm_co2,
1614     $ zzlev,zzlay,
1615     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
1616     $ zplay,zplev,pphi,zpopsk,
1617     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
1618     $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
1619
1620         DO l=1,nlayer
1621           DO ig=1,ngrid
1622              pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
1623              pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
1624              pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
1625              q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
1626           ENDDO
1627        ENDDO
1628
1629        DO ig=1,ngrid
1630          q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
1631        ENDDO
1632
1633        DO iq=1,nq
1634         DO l=1,nlayer
1635           DO ig=1,ngrid
1636             pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
1637           ENDDO
1638         ENDDO
1639        ENDDO
1640
1641        lmax_th_out(:)=real(lmax_th(:))
1642
1643      else   !of if calltherm
1644        lmax_th(:)=0
1645        wstar(:)=0.
1646        hfmax_th(:)=0.
1647        lmax_th_out(:)=0.
1648      end if
1649
1650c-----------------------------------------------------------------------
1651c   7. Dry convective adjustment:
1652c   -----------------------------
1653
1654      IF(calladj) THEN
1655
1656         DO l=1,nlayer
1657            DO ig=1,ngrid
1658               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
1659            ENDDO
1660         ENDDO
1661         zduadj(:,:)=0
1662         zdvadj(:,:)=0
1663         zdhadj(:,:)=0
1664
1665         CALL convadj(ngrid,nlayer,nq,ptimestep,
1666     $                zplay,zplev,zpopsk,lmax_th,
1667     $                pu,pv,zh,pq,
1668     $                pdu,pdv,zdh,pdq,
1669     $                zduadj,zdvadj,zdhadj,
1670     $                zdqadj)
1671
1672         DO l=1,nlayer
1673            DO ig=1,ngrid
1674               pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
1675               pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
1676               pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
1677
1678               zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1679            ENDDO
1680         ENDDO
1681
1682           DO iq=1, nq
1683            DO l=1,nlayer
1684              DO ig=1,ngrid
1685                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
1686              ENDDO
1687            ENDDO
1688           ENDDO
1689      ENDIF ! of IF(calladj)
1690
1691c-----------------------------------------------------
1692c    8. Non orographic Gravity waves :
1693c -------------------------------------------------
1694
1695      IF (calllott_nonoro) THEN
1696
1697         CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,
1698     &               cpnew,rnew,
1699     &               zplay,
1700     &               zmax_th,                      ! max altitude reached by thermals (m)
1701     &               pt, pu, pv,
1702     &               pdt, pdu, pdv,
1703     &               zustrhi,zvstrhi,
1704     &               d_t_hin, d_u_hin, d_v_hin)
1705         IF (calljliu_gwimix) THEN
1706          CALL nonoro_gwd_mix(ngrid,nlayer,ptimestep,
1707     &               nq,cpnew, rnew,
1708     &               zplay,
1709     &               zmax_th,
1710     &               pt, pu, pv, pq, zh,
1711                !loss,  chemical reaction loss rates
1712     &               pdt, pdu, pdv, pdq, zdh,
1713                !  zustrhi,zvstrhi,
1714     &               zdq_mix, d_t_mix, d_u_mix, d_v_mix)
1715         ENDIF
1716
1717!  Update tendencies
1718         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)
1719     &                         +d_t_hin(1:ngrid,1:nlayer)
1720         pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
1721     &                         +d_u_hin(1:ngrid,1:nlayer)
1722         pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
1723     &                         +d_v_hin(1:ngrid,1:nlayer)
1724!  Update tendencies of gw mixing
1725        IF (calljliu_gwimix) THEN
1726         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)
1727     &                         +d_t_mix(1:ngrid,1:nlayer)
1728         pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
1729     &                         +d_u_mix(1:ngrid,1:nlayer)
1730         pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
1731     &                         +d_v_mix(1:ngrid,1:nlayer)
1732       pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)
1733     &                       +zdq_mix(1:ngrid,1:nlayer,1:nq)
1734        ENDIF
1735
1736
1737      ENDIF ! of IF (calllott_nonoro)
1738
1739c-----------------------------------------------------------------------
1740c   9. Specific parameterizations for tracers
1741c:   -----------------------------------------
1742
1743
1744c   9a. Water and ice
1745c     ---------------
1746
1747c        ---------------------------------------
1748c        Water ice condensation in the atmosphere
1749c        ----------------------------------------
1750         IF (water) THEN
1751
1752           call watercloud(ngrid,nlayer,ptimestep,
1753     &                zplev,zplay,pdpsrf,zzlay, pt,pdt,
1754     &                pq,pdq,zdqcloud,zdtcloud,
1755     &                nq,tau,tauscaling,rdust,rice,nuice,
1756     &                rsedcloud,rhocloud,totcloudfrac)
1757c Temperature variation due to latent heat release
1758           if (activice) then
1759               pdt(1:ngrid,1:nlayer) =
1760     &          pdt(1:ngrid,1:nlayer) +
1761     &          zdtcloud(1:ngrid,1:nlayer)
1762           endif
1763
1764! increment water vapour and ice atmospheric tracers tendencies
1765           pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
1766     &       pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
1767     &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
1768           pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1769     &       pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1770     &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
1771
1772           if (hdo) then
1773! increment HDO vapour and ice atmospheric tracers tendencies
1774           pdq(1:ngrid,1:nlayer,igcm_hdo_vap) =
1775     &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) +
1776     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
1777           pdq(1:ngrid,1:nlayer,igcm_hdo_ice) =
1778     &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) +
1779     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
1780           endif !hdo
1781
1782! increment dust and ccn masses and numbers
1783! We need to check that we have Nccn & Ndust > 0
1784! This is due to single precision rounding problems
1785           if (microphys) then
1786             pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1787     &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1788     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
1789             pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1790     &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1791     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
1792             where (pq(:,:,igcm_ccn_mass) +
1793     &       ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1794               pdq(:,:,igcm_ccn_mass) =
1795     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1796               pdq(:,:,igcm_ccn_number) =
1797     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1798             end where
1799             where (pq(:,:,igcm_ccn_number) +
1800     &       ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1801               pdq(:,:,igcm_ccn_mass) =
1802     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1803               pdq(:,:,igcm_ccn_number) =
1804     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1805             end where
1806           endif
1807
1808           if (scavenging) then
1809             pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1810     &         pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1811     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
1812             pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1813     &         pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1814     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
1815             where (pq(:,:,igcm_dust_mass) +
1816     &       ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1817               pdq(:,:,igcm_dust_mass) =
1818     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1819               pdq(:,:,igcm_dust_number) =
1820     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1821             end where
1822             where (pq(:,:,igcm_dust_number) +
1823     &       ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1824               pdq(:,:,igcm_dust_mass) =
1825     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1826               pdq(:,:,igcm_dust_number) =
1827     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1828             end where
1829           endif ! of if scavenging
1830
1831         END IF  ! of IF (water)
1832
1833c   9a bis. CO2 clouds (CL & JA)
1834c        ---------------------------------------
1835c        CO2 ice cloud condensation in the atmosphere
1836c        ----------------------------------------
1837c flag needed in callphys.def:
1838c               co2clouds=.true. is mandatory (default is .false.)
1839c               co2useh2o=.true. if you want to allow co2 condensation
1840c                                on water ice particles
1841c               meteo_flux=.true. if you want to add a meteoritic
1842c                                 supply of CCN
1843c               CLFvaryingCO2=.true. if you want to have a sub-grid
1844c                                    temperature distribution
1845c               spantCO2=integer (i.e. 3) amplitude of the sub-grid T disti
1846c               nuiceco2_sed=0.2 variance of the size distribution for the
1847c                                sedimentation
1848c               nuiceco2_ref=0.2 variance of the size distribution for the
1849c                                nucleation
1850c               imicroco2=50     micro-timestep is 1/50 of physical timestep
1851        zdqssed_co2(:) = 0.
1852        zdqssed_ccn(:,:) = 0.
1853
1854         IF (co2clouds) THEN
1855           call co2cloud(ngrid,nlayer,ptimestep,
1856     &           zplev,zplay,pdpsrf,zzlay,pt,pdt,
1857     &           pq,pdq,zdqcloudco2,zdtcloudco2,
1858     &           nq,tau,tauscaling,rdust,rice,riceco2,nuice,
1859     &           rhocloud, rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2,
1860     &           zdqssed_ccn,pdu,pu,zcondicea_co2microp)
1861
1862          DO iq=1, nq
1863            DO ig=1,ngrid
1864             DO islope = 1,nslope
1865              dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+
1866     &         zdqssed_ccn(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
1867             ENDDO !(islope)
1868            ENDDO  ! (ig)
1869           ENDDO    ! (iq)q)
1870c Temperature variation due to latent heat release
1871               pdt(1:ngrid,1:nlayer) =
1872     &              pdt(1:ngrid,1:nlayer) +
1873     &              zdtcloudco2(1:ngrid,1:nlayer)
1874
1875! increment dust and ccn masses and numbers
1876! We need to check that we have Nccn & Ndust > 0
1877! This is due to single precision rounding problems
1878! increment dust tracers tendancies
1879               pdq(:,:,igcm_dust_mass) = pdq(:,:,igcm_dust_mass)
1880     &                                 + zdqcloudco2(:,:,igcm_dust_mass)
1881
1882               pdq(:,:,igcm_dust_number) = pdq(:,:,igcm_dust_number)
1883     &                               + zdqcloudco2(:,:,igcm_dust_number)
1884
1885               pdq(:,:,igcm_co2) = pdq(:,:,igcm_co2)
1886     &                             + zdqcloudco2(:,:,igcm_co2)
1887
1888               pdq(:,:,igcm_co2_ice) = pdq(:,:,igcm_co2_ice)
1889     &                                 + zdqcloudco2(:,:,igcm_co2_ice)
1890
1891               pdq(:,:,igcm_ccnco2_mass) = pdq(:,:,igcm_ccnco2_mass)
1892     &                               + zdqcloudco2(:,:,igcm_ccnco2_mass)
1893
1894               pdq(:,:,igcm_ccnco2_number) = pdq(:,:,igcm_ccnco2_number)
1895     &                             + zdqcloudco2(:,:,igcm_ccnco2_number)
1896
1897             if (meteo_flux) then
1898               pdq(:,:,igcm_ccnco2_meteor_mass) =
1899     &                 pdq(:,:,igcm_ccnco2_meteor_mass) +
1900     &                 zdqcloudco2(:,:,igcm_ccnco2_meteor_mass)
1901
1902               pdq(:,:,igcm_ccnco2_meteor_number) =
1903     &                 pdq(:,:,igcm_ccnco2_meteor_number)
1904     &                 + zdqcloudco2(:,:,igcm_ccnco2_meteor_number)
1905             end if
1906!Update water ice clouds values as well
1907             if (co2useh2o) then
1908                pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1909     &               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1910     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
1911                pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1912     &               pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1913     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
1914                pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1915     &               pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1916     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
1917
1918               pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1919     &                pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
1920     &                zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice)
1921
1922               pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1923     &                pdq(:,:,igcm_ccnco2_h2o_mass_ccn) +
1924     &                zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn)
1925
1926               pdq(:,:,igcm_ccnco2_h2o_number) =
1927     &                   pdq(:,:,igcm_ccnco2_h2o_number) +
1928     &                   zdqcloudco2(:,:,igcm_ccnco2_h2o_number)
1929
1930c     Negative values?
1931                where (pq(:,:,igcm_ccn_mass) +
1932     &                 ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1933                  pdq(:,:,igcm_ccn_mass) =
1934     &               - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1935                  pdq(:,:,igcm_ccn_number) =
1936     &               - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1937                end where
1938c     Negative values?
1939                where (pq(:,:,igcm_ccn_number) +
1940     &                 ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1941                  pdq(:,:,igcm_ccn_mass) =
1942     &              - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1943                  pdq(:,:,igcm_ccn_number) =
1944     &              - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1945                end where
1946             where (pq(:,:,igcm_ccnco2_h2o_mass_ice) +
1947     &              pq(:,:,igcm_ccnco2_h2o_mass_ccn) +
1948     &              (pdq(:,:,igcm_ccnco2_h2o_mass_ice) +
1949     &               pdq(:,:,igcm_ccnco2_h2o_mass_ccn)
1950     &               )*ptimestep < 0.)
1951               pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1952     &            - pq(:,:,igcm_ccnco2_h2o_mass_ice)
1953     &               /ptimestep + 1.e-30
1954               pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1955     &            - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
1956     &               /ptimestep + 1.e-30
1957               pdq(:,:,igcm_ccnco2_h2o_number) =
1958     &            - pq(:,:,igcm_ccnco2_h2o_number)
1959     &                /ptimestep + 1.e-30
1960             end where
1961
1962             where (pq(:,:,igcm_ccnco2_h2o_number) +
1963     &              (pdq(:,:,igcm_ccnco2_h2o_number)
1964     &               )*ptimestep < 0.)
1965               pdq(:,:,igcm_ccnco2_h2o_mass_ice) =
1966     &            - pq(:,:,igcm_ccnco2_h2o_mass_ice)
1967     &               /ptimestep + 1.e-30
1968               pdq(:,:,igcm_ccnco2_h2o_mass_ccn) =
1969     &            - pq(:,:,igcm_ccnco2_h2o_mass_ccn)
1970     &               /ptimestep + 1.e-30
1971               pdq(:,:,igcm_ccnco2_h2o_number) =
1972     &            - pq(:,:,igcm_ccnco2_h2o_number)
1973     &                /ptimestep + 1.e-30
1974             end where
1975             endif ! of if (co2useh2o)
1976c     Negative values?
1977             where (pq(:,:,igcm_ccnco2_mass) +
1978     &              ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
1979               pdq(:,:,igcm_ccnco2_mass) =
1980     &            - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1981               pdq(:,:,igcm_ccnco2_number) =
1982     &            - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1983             end where
1984             where (pq(:,:,igcm_ccnco2_number) +
1985     &              ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
1986               pdq(:,:,igcm_ccnco2_mass) =
1987     &            - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
1988               pdq(:,:,igcm_ccnco2_number) =
1989     &            - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
1990             end where
1991
1992c     Negative values?
1993             where (pq(:,:,igcm_dust_mass) +
1994     &              ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1995               pdq(:,:,igcm_dust_mass) =
1996     &           - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1997               pdq(:,:,igcm_dust_number) =
1998     &           - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1999             end where
2000             where (pq(:,:,igcm_dust_number) +
2001     &              ptimestep*pdq(:,:,igcm_dust_number) < 0.)
2002               pdq(:,:,igcm_dust_mass) =
2003     &           - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
2004               pdq(:,:,igcm_dust_number) =
2005     &           - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
2006             end where
2007             if (meteo_flux) then
2008             where (pq(:,:,igcm_ccnco2_meteor_mass) +
2009     &              ptimestep*pdq(:,:,igcm_ccnco2_meteor_mass) < 0.)
2010               pdq(:,:,igcm_ccnco2_meteor_mass) =
2011     &            - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30
2012               pdq(:,:,igcm_ccnco2_meteor_number) =
2013     &            - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30
2014             end where
2015             where (pq(:,:,igcm_ccnco2_meteor_number) +
2016     &              ptimestep*pdq(:,:,igcm_ccnco2_meteor_number) < 0.)
2017               pdq(:,:,igcm_ccnco2_meteor_mass) =
2018     &            - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30
2019               pdq(:,:,igcm_ccnco2_meteor_number) =
2020     &            - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30
2021             end where
2022             end if
2023      END IF ! of IF (co2clouds)
2024
2025c   9b. Aerosol particles
2026c     -------------------
2027c        ----------
2028c        Dust devil :
2029c        ----------
2030         IF(callddevil) then
2031           call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
2032     &                zdqdev,zdqsdev)
2033
2034           if (dustbin.ge.1) then
2035              do iq=1,nq
2036                 DO l=1,nlayer
2037                    DO ig=1,ngrid
2038                       pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
2039                    ENDDO
2040                 ENDDO
2041              enddo
2042              do iq=1,nq
2043                 DO ig=1,ngrid
2044                   DO islope = 1,nslope
2045                    dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
2046     &          zdqsdev(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
2047                   ENDDO
2048                 ENDDO
2049              enddo
2050           endif  ! of if (dustbin.ge.1)
2051
2052         END IF ! of IF (callddevil)
2053
2054c        -------------
2055c        Sedimentation :   acts also on water ice
2056c        -------------
2057         IF (sedimentation) THEN
2058           zdqsed(1:ngrid,1:nlayer,1:nq)=0
2059           zdqssed(1:ngrid,1:nq)=0
2060
2061c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
2062c Zdqssed isn't
2063
2064           call callsedim(ngrid,nlayer,ptimestep,
2065     &                zplev,zzlev,zzlay,pt,pdt,
2066     &                rdust,rstormdust,rtopdust,
2067     &                rice,rsedcloud,rhocloud,
2068     &                pq,pdq,zdqsed,zdqssed,nq,
2069     &                tau,tauscaling)
2070
2071
2072c Flux at the surface of co2 ice computed in co2cloud microtimestep
2073           IF (rdstorm) THEN
2074c             Storm dust cannot sediment to the surface
2075              DO ig=1,ngrid
2076                 zdqsed(ig,1,igcm_stormdust_mass)=
2077     &             zdqsed(ig,1,igcm_stormdust_mass)+
2078     &             zdqssed(ig,igcm_stormdust_mass) /
2079     &             ((pplev(ig,1)-pplev(ig,2))/g)
2080                 zdqsed(ig,1,igcm_stormdust_number)=
2081     &             zdqsed(ig,1,igcm_stormdust_number)+
2082     &             zdqssed(ig,igcm_stormdust_number) /
2083     &               ((pplev(ig,1)-pplev(ig,2))/g)
2084                 zdqssed(ig,igcm_stormdust_mass)=0.
2085                 zdqssed(ig,igcm_stormdust_number)=0.
2086              ENDDO
2087           ENDIF !rdstorm
2088
2089           DO iq=1, nq
2090             DO l=1,nlayer
2091               DO ig=1,ngrid
2092                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
2093               ENDDO
2094             ENDDO
2095           ENDDO
2096           DO iq=1, nq
2097             DO ig=1,ngrid
2098               DO islope = 1,nslope
2099                    dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
2100     &          zdqssed(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
2101               ENDDO
2102             ENDDO
2103           ENDDO
2104
2105         END IF   ! of IF (sedimentation)
2106
2107c Add lifted dust to tendancies after sedimentation in the LES (AC)
2108      IF (turb_resolved) THEN
2109           DO iq=1, nq
2110            DO l=1,nlayer
2111              DO ig=1,ngrid
2112                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
2113              ENDDO
2114            ENDDO
2115           ENDDO
2116           DO iq=1, nq
2117              DO ig=1,ngrid
2118                 dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
2119              ENDDO
2120           ENDDO
2121      ENDIF
2122c
2123c   9c. Chemical species
2124c     ------------------
2125
2126#ifndef MESOSCALE
2127c        --------------
2128c        photochemistry :
2129c        --------------
2130         IF (photochem) then
2131
2132           if (modulo(icount-1,ichemistry).eq.0) then
2133           ! compute chemistry every ichemistry physics step
2134
2135!           dust and ice surface area
2136            call surfacearea(ngrid, nlayer, naerkind,
2137     $                       ptimestep, zplay, zzlay,
2138     $                       pt, pq, pdq, nq,
2139     $                       rdust, rice, tau, tauscaling,
2140     $                       surfdust, surfice)
2141!           call photochemistry
2142            DO ig = 1,ngrid
2143              qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig))
2144            ENDDO
2145            call calchim(ngrid,nlayer,nq,
2146     &                   ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
2147     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
2148     $                   zdqcloud,zdqscloud,tau(:,1),
2149     $                   qsurf_tmp(:,igcm_co2),
2150     $                   pu,pdu,pv,pdv,surfdust,surfice)
2151
2152            endif ! of if (modulo(icount-1,ichemistry).eq.0)
2153
2154           ! increment values of tracers:
2155           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
2156                      ! tracers is zero anyways
2157             DO l=1,nlayer
2158               DO ig=1,ngrid
2159                 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
2160               ENDDO
2161             ENDDO
2162           ENDDO ! of DO iq=1,nq
2163
2164           ! add condensation tendency for H2O2
2165           if (igcm_h2o2.ne.0) then
2166             DO l=1,nlayer
2167               DO ig=1,ngrid
2168                 pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
2169     &                                +zdqcloud(ig,l,igcm_h2o2)
2170               ENDDO
2171             ENDDO
2172           endif
2173
2174           ! increment surface values of tracers:
2175           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
2176                      ! tracers is zero anyways
2177             DO ig=1,ngrid
2178              DO islope = 1,nslope
2179               dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) +
2180     &           zdqschim(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
2181              ENDDO
2182             ENDDO
2183           ENDDO ! of DO iq=1,nq
2184
2185           ! add condensation tendency for H2O2
2186           if (igcm_h2o2.ne.0) then
2187             DO ig=1,ngrid
2188              DO islope = 1,nslope
2189              dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+
2190     &      zdqscloud(ig,igcm_h2o2)*cos(pi*def_slope_mean(islope)/180.)
2191              ENDDO
2192             ENDDO
2193           endif
2194
2195         END IF  ! of IF (photochem)
2196#endif
2197
2198
2199#ifndef MESOSCALE
2200c-----------------------------------------------------------------------
2201c   10. THERMOSPHERE CALCULATION
2202c-----------------------------------------------------------------------
2203
2204      if (callthermos) then
2205        call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
2206     $     mu0,ptimestep,ptime,zday,tsurf_meshavg,zzlev,zzlay,
2207     &     pt,pq,pu,pv,pdt,pdq,
2208     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
2209     $     PhiEscH,PhiEscH2,PhiEscD)
2210
2211        DO l=1,nlayer
2212          DO ig=1,ngrid
2213            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
2214            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l)
2215            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
2216            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
2217            DO iq=1, nq
2218              pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
2219            ENDDO
2220          ENDDO
2221        ENDDO
2222
2223      endif ! of if (callthermos)
2224#endif
2225
2226c-----------------------------------------------------------------------
2227c   11. Carbon dioxide condensation-sublimation:
2228c     (should be the last atmospherical physical process to be computed)
2229c   -------------------------------------------
2230      IF (tituscap) THEN
2231         !!! get the actual co2 seasonal cap from Titus observations
2232         CALL geticecover(ngrid, 180.*zls/pi,
2233     .                  180.*longitude/pi, 180.*latitude/pi,
2234     .                  qsurf_tmp(:,igcm_co2) )
2235         qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000.
2236      ENDIF
2237
2238
2239      IF (callcond) THEN
2240         zdtc(:,:) = 0.
2241         zdtsurfc(:,:) = 0.
2242         zduc(:,:) = 0.
2243         zdvc(:,:) = 0.
2244         zdqc(:,:,:) = 0.
2245         zdqsc(:,:,:) = 0.
2246         CALL co2condens(ngrid,nlayer,nq,nslope,ptimestep,
2247     $              capcal,zplay,zplev,tsurf,pt,
2248     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
2249     $              qsurf(:,igcm_co2,:),perennial_co2ice,
2250     $              albedo,emis,rdust,
2251     $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
2252     $              fluxsurf_dn_sw,zls,
2253     $              zdqssed_co2,zcondicea_co2microp,
2254     &              zdqsc)
2255
2256          if (ngrid == 1) then ! For the 1D model
2257          ! CO2cond_ps is a coefficient to control the surface pressure change
2258              pdpsrf = CO2cond_ps*pdpsrf
2259              zduc   = CO2cond_ps*zduc
2260              zdvc   = CO2cond_ps*zdvc
2261              zdqc   = CO2cond_ps*zdqc
2262          endif
2263
2264         DO iq=1, nq
2265           DO ig=1,ngrid
2266              dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:)
2267           ENDDO  ! (ig)
2268         ENDDO    ! (iq)
2269         DO l=1,nlayer
2270           DO ig=1,ngrid
2271             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
2272             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
2273             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
2274           ENDDO
2275         ENDDO
2276         DO ig=1,ngrid
2277           zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:)
2278         ENDDO
2279
2280           DO iq=1, nq
2281            DO l=1,nlayer
2282              DO ig=1,ngrid
2283                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
2284              ENDDO
2285            ENDDO
2286           ENDDO
2287
2288#ifndef MESOSCALE
2289        ! update surface pressure
2290        DO ig=1,ngrid
2291          ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
2292        ENDDO
2293        ! update pressure levels
2294        DO l=1,nlayer
2295         DO ig=1,ngrid
2296          zplay(ig,l) = aps(l) + bps(l)*ps(ig)
2297          zplev(ig,l) = ap(l)  + bp(l)*ps(ig)
2298         ENDDO
2299        ENDDO
2300        zplev(:,nlayer+1) = 0.
2301
2302    ! Calculation of zzlay and zzlay with udpated pressure and temperature
2303        DO ig=1,ngrid
2304         zzlay(ig,1)=-(log(zplay(ig,1)/ps(ig)))*rnew(ig,1)*
2305     &      (pt(ig,1)+pdt(ig,1)*ptimestep) /g
2306
2307          DO l=2,nlayer
2308
2309           ! compute "mean" temperature of the layer
2310            if((pt(ig,l)+pdt(ig,l)*ptimestep) .eq.
2311     &               (pt(ig,l-1)+pdt(ig,l-1)*ptimestep)) then
2312               tlaymean= pt(ig,l)+pdt(ig,l)*ptimestep
2313            else
2314               tlaymean=((pt(ig,l)+pdt(ig,l)*ptimestep)-
2315     &        (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))/
2316     &        log((pt(ig,l)+pdt(ig,l)*ptimestep)/
2317     &        (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))
2318            endif
2319
2320            ! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
2321            gz(ig,l)=g*(rad**2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))**2
2322
2323
2324            zzlay(ig,l)=zzlay(ig,l-1)-
2325     &     (log(zplay(ig,l)/zplay(ig,l-1))*rnew(ig,l)*tlaymean/gz(ig,l))
2326
2327
2328        !   update layers altitude
2329             z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
2330             z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
2331             zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
2332          ENDDO
2333        ENDDO
2334#endif
2335      ENDIF  ! of IF (callcond)
2336
2337c-----------------------------------------------------------------------
2338c  Updating tracer budget on surface
2339c-----------------------------------------------------------------------
2340        DO iq=1, nq
2341          DO ig=1,ngrid
2342           DO islope = 1,nslope
2343            qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+
2344     &                ptimestep*dqsurf(ig,iq,islope)
2345           ENDDO
2346          ENDDO  ! (ig)
2347        ENDDO    ! (iq)
2348c-----------------------------------------------------------------------
2349c   12. Surface  and sub-surface soil temperature
2350c-----------------------------------------------------------------------
2351c
2352c
2353c   12.1 Increment Surface temperature:
2354c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2355
2356      DO ig=1,ngrid
2357        DO islope = 1,nslope
2358            tsurf(ig,islope)=tsurf(ig,islope)+
2359     &            ptimestep*zdtsurf(ig,islope)
2360       ENDDO
2361      ENDDO
2362
2363c  Prescribe a cold trap at south pole (except at high obliquity !!)
2364c  Temperature at the surface is set there to be the temperature
2365c  corresponding to equilibrium temperature between phases of CO2
2366
2367
2368      IF (water) THEN
2369!#ifndef MESOSCALE
2370!         if (caps.and.(obliquit.lt.27.)) then => now done in co2condens
2371           ! NB: Updated surface pressure, at grid point 'ngrid', is
2372           !     ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
2373!           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
2374!     &                     (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
2375!           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*ps(ngrid)))
2376!         endif
2377!#endif
2378c       -------------------------------------------------------------
2379c       Change of surface albedo in case of ground frost
2380c       everywhere except on the north permanent cap and in regions
2381c       covered by dry ice.
2382c              ALWAYS PLACE these lines after co2condens !!!
2383c       -------------------------------------------------------------
2384         do ig = 1,ngrid
2385          do islope = 1,nslope
2386           if (abs(qsurf(ig,igcm_co2,islope)) < 1.e-10) then ! No CO2 frost
2387
2388             if (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)
2389     & then ! There is H2O frost
2390               if (cst_cap_albedo .and. watercaptag(ig) .and.
2391     & abs(perennial_co2ice(ig,islope)) < 1.e-10) then ! Water cap remains unchanged by water frost deposition and no CO2 perennial ice
2392                 albedo(ig,:,islope) = albedo_h2o_cap
2393                 emis(ig,islope) = 1.
2394               else
2395                 albedo(ig,:,islope) = albedo_h2o_frost
2396                 emis(ig,islope) = 1.
2397               endif
2398             else ! No H2O frost
2399               if (abs(perennial_co2ice(ig,islope)) < 1.e-10 .and.
2400     & watercaptag(ig)) then ! No CO2 perennial ice but there is water cap
2401                 albedo(ig,:,islope) = albedo_h2o_cap
2402                 emis(ig,islope) = 1.
2403               endif
2404             endif
2405
2406           endif
2407          enddo ! islope
2408         enddo  ! of do ig=1,ngrid
2409      ENDIF  ! of IF (water)
2410
2411c
2412c   12.2 Compute soil temperatures and subsurface heat flux:
2413c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2414      IF (callsoil) THEN
2415c       Thermal inertia feedback
2416        IF (surfaceice_tifeedback.or.poreice_tifeedback) THEN
2417       
2418           CALL waterice_tifeedback(ngrid,nsoilmx,nslope,
2419     s               qsurf(:,igcm_h2o_ice,:),pore_icefraction,
2420     s               inertiesoil_tifeedback(:,:,:))
2421     
2422         CALL soil(ngrid,nsoilmx,.false.,inertiesoil_tifeedback,
2423     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
2424        ELSE
2425         CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
2426     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
2427        ENDIF
2428      ENDIF
2429
2430c     To avoid negative values
2431      IF (rdstorm) THEN
2432           where (pq(:,:,igcm_stormdust_mass) +
2433     &      ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.)
2434             pdq(:,:,igcm_stormdust_mass) =
2435     &       - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
2436             pdq(:,:,igcm_stormdust_number) =
2437     &       - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30
2438           end where
2439           where (pq(:,:,igcm_stormdust_number) +
2440     &      ptimestep*pdq(:,:,igcm_stormdust_number) < 0.)
2441             pdq(:,:,igcm_stormdust_mass) =
2442     &       - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
2443             pdq(:,:,igcm_stormdust_number) =
2444     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
2445           end where
2446
2447            where (pq(:,:,igcm_dust_mass) +
2448     &      ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
2449             pdq(:,:,igcm_dust_mass) =
2450     &       - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
2451             pdq(:,:,igcm_dust_number) =
2452     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
2453           end where
2454           where (pq(:,:,igcm_dust_number) +
2455     &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
2456             pdq(:,:,igcm_dust_mass) =
2457     &       - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
2458             pdq(:,:,igcm_dust_number) =
2459     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
2460           end where
2461      ENDIF !(rdstorm)
2462
2463c-----------------------------------------------------------------------
2464c   J. Naar : Surface and sub-surface water ice
2465c-----------------------------------------------------------------------
2466c
2467c
2468c   Increment Watercap (surface h2o reservoirs):
2469c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2470
2471      DO ig=1,ngrid
2472        DO islope = 1,nslope
2473         watercap(ig,islope)=watercap(ig,islope)+
2474     s   ptimestep*dwatercap(ig,islope)
2475        ENDDO
2476      ENDDO
2477
2478      IF (refill_watercap) THEN
2479
2480        DO ig = 1,ngrid
2481         DO islope = 1,nslope
2482           if (watercaptag(ig) .and. (qsurf(ig,igcm_h2o_ice,islope)
2483     &        > frost_metam_threshold)) then
2484
2485                watercap(ig,islope) = watercap(ig,islope)
2486     &          + qsurf(ig,igcm_h2o_ice,islope)
2487     &          - frost_metam_threshold
2488                qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold
2489           endif ! watercaptag
2490          ENDDO
2491        ENDDO
2492
2493      ENDIF ! refill_watercap
2494
2495c-----------------------------------------------------------------------
2496c  13. Write output files
2497c  ----------------------
2498      call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
2499     &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
2500
2501c    -------------------------------
2502c    Dynamical fields incrementation
2503c    -------------------------------
2504c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
2505      ! temperature, zonal and meridional wind
2506      DO l=1,nlayer
2507        DO ig=1,ngrid
2508          zt(ig,l)=pt(ig,l)  + pdt(ig,l)*ptimestep
2509          zu(ig,l)=pu(ig,l)  + pdu(ig,l)*ptimestep
2510          zv(ig,l)=pv(ig,l)  + pdv(ig,l)*ptimestep
2511        ENDDO
2512      ENDDO
2513
2514      ! tracers
2515      DO iq=1, nq
2516        DO l=1,nlayer
2517          DO ig=1,ngrid
2518            zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
2519          ENDDO
2520        ENDDO
2521      ENDDO
2522
2523      ! Density
2524      DO l=1,nlayer
2525         DO ig=1,ngrid
2526            rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
2527         ENDDO
2528      ENDDO
2529
2530      ! Potential Temperature
2531
2532       DO ig=1,ngrid
2533          DO l=1,nlayer
2534              zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp
2535          ENDDO
2536       ENDDO
2537
2538c    Compute surface stress : (NB: z0 is a common in surfdat.h)
2539c     DO ig=1,ngrid
2540c        cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
2541c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
2542c     ENDDO
2543
2544c     Sum of fluxes in solar spectral bands (for output only)
2545      fluxtop_dn_sw_tot(1:ngrid)=fluxtop_dn_sw(1:ngrid,1) +
2546     &                           fluxtop_dn_sw(1:ngrid,2)
2547      fluxtop_up_sw_tot(1:ngrid)=fluxtop_up_sw(1:ngrid,1) +
2548     &                           fluxtop_up_sw(1:ngrid,2)
2549      fluxsurf_dn_sw_tot(1:ngrid,1:nslope)=
2550     &                           fluxsurf_dn_sw(1:ngrid,1,1:nslope) +
2551     &                           fluxsurf_dn_sw(1:ngrid,2,1:nslope)
2552      fluxsurf_up_sw_tot(1:ngrid)=fluxsurf_up_sw(1:ngrid,1) +
2553     &                            fluxsurf_up_sw(1:ngrid,2)
2554
2555c ******* TEST ******************************************************
2556      ztim1 = 999
2557      DO l=1,nlayer
2558        DO ig=1,ngrid
2559           if (pt(ig,l).lt.ztim1) then
2560               ztim1 = pt(ig,l)
2561               igmin = ig
2562               lmin = l
2563           end if
2564        ENDDO
2565      ENDDO
2566      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
2567        write(*,*) 'PHYSIQ: stability WARNING :'
2568        write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
2569     &              'ig l =', igmin, lmin
2570      end if
2571
2572c        ----------------------------------------------------------
2573c        ----------------------------------------------------------
2574c        INTERPOLATIONS IN THE SURFACE-LAYER
2575c        ----------------------------------------------------------
2576c        ----------------------------------------------------------
2577
2578           n_out=0 ! number of elements in the z_out array.
2579                   ! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
2580                   ! to 5
2581           IF (n_out .ne. 0) THEN
2582
2583           IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
2584           IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
2585           IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
2586
2587           z_out(:)=[3.,2.,1.,0.5,0.1]
2588           u_out(:,:)=0.
2589           T_out(:,:)=0.
2590
2591           call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
2592     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,q2,tsurf(:,iflat),
2593     & zh,zq(:,1,igcm_h2o_vap),qsurf(:,igcm_h2o_ice,iflat),mmean(:,1),
2594     & z_out,n_out,T_out,u_out,ustar,tstar,vhf,vvv)
2595                   ! pourquoi ustar recalcule ici? fait dans vdifc.
2596
2597#ifndef MESOSCALE
2598            DO n=1,n_out
2599               write(zstring, '(F8.6)') z_out(n)
2600               call write_output('T_out_'//trim(zstring),
2601     &   'potential temperature at z_out','K',T_out(:,n))
2602               call write_output('u_out_'//trim(zstring),
2603     &   'horizontal velocity norm at z_out','m/s',u_out(:,n))
2604            ENDDO
2605            call write_output('u_star',
2606     &   'friction velocity','m/s',ustar)
2607            call write_output('teta_star',
2608     &   'friction potential temperature','K',tstar)
2609            call write_output('vvv',
2610     &   'Vertical velocity variance at zout','m',vvv)
2611            call write_output('vhf',
2612     &   'Vertical heat flux at zout','m',vhf)
2613#else
2614         T_out1(:)=T_out(:,1)
2615         u_out1(:)=u_out(:,1)
2616#endif
2617
2618         ENDIF
2619
2620c        ----------------------------------------------------------
2621c        ----------------------------------------------------------
2622c        END OF SURFACE LAYER INTERPOLATIONS
2623c        ----------------------------------------------------------
2624c        ----------------------------------------------------------
2625
2626#ifndef MESOSCALE
2627c        -------------------------------------------------------------------
2628c        Writing NetCDF file  "RESTARTFI" at the end of the run
2629c        -------------------------------------------------------------------
2630c        Note: 'restartfi' is stored just before dynamics are stored
2631c              in 'restart'. Between now and the writting of 'restart',
2632c              there will have been the itau=itau+1 instruction and
2633c              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
2634c              thus we store for time=time+dtvr
2635
2636         ! default: not writing a restart file at this time step
2637         write_restart=.false.
2638         IF (ecritstart.GT.0) THEN
2639           ! For when we store multiple time steps in the restart file
2640           IF (MODULO(icount*iphysiq,ecritstart).EQ.0) THEN
2641             write_restart=.true.
2642           ENDIF
2643         ENDIF
2644         IF (lastcall) THEN
2645           ! Always write a restart at the end of the simulation
2646           write_restart=.true.
2647         ENDIF
2648
2649          IF (write_restart) THEN
2650           IF (grid_type==unstructured) THEN !IF DYNAMICO
2651
2652             ! When running Dynamico, no need to add a dynamics time step to ztime_fin
2653             IF (ptime.LE. 1.E-10) THEN
2654             ! Residual ptime occurs with Dynamico
2655             ztime_fin = pday !+ ptime + ptimestep/(float(iphysiq)*daysec)
2656     .               - day_ini - time_phys
2657             ELSE
2658             ztime_fin = pday + ptime  !+ ptimestep/(float(iphysiq)*daysec)
2659     .                  - day_ini - time_phys
2660             ENDIF
2661             if (ecritstart==0) then
2662                ztime_fin = ztime_fin-(day_end-day_ini)
2663             endif
2664
2665           ELSE ! IF LMDZ
2666
2667            if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change
2668              ztime_fin = pday - day_ini + ptime
2669     &                    + ptimestep/(float(iphysiq)*daysec)
2670            else !IF ONE RESTART final time in top of day_end
2671          ztime_fin = pday - day_ini-(day_end-day_ini)
2672     &                    + ptime  + ptimestep/(float(iphysiq)*daysec)
2673            endif
2674
2675           ENDIF ! of IF (grid_type==unstructured)
2676          write(*,'(A,I7,A,F12.5)')
2677     .         'PHYSIQ: writing in restartfi ; icount=',
2678     .          icount,' date=',ztime_fin
2679
2680          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,
2681     .                ptimestep,ztime_fin,
2682     .                tsurf,tsoil,inertiesoil,albedo,
2683     .                emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,
2684     .                wstar,watercap,perennial_co2ice)
2685          ENDIF ! of IF (write_restart)
2686
2687#endif
2688
2689c         IF (ngrid.NE.1) then
2690
2691c        -------------------------------------------------------------------
2692c        Calculation of diagnostic variables written in both stats and
2693c          diagfi files
2694c        -------------------------------------------------------------------
2695         do ig=1,ngrid
2696       if(mu0(ig).le.0.01) then
2697        fluxsurf_dir_dn_sw(ig) = 0.
2698       else
2699            if (water) then
2700             ! both water and dust contribute
2701            fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_sol*mu0(ig)*
2702     &                    exp(-(tau(ig,iaer_dust_doubleq)+
2703     &                          tau(ig,iaer_h2o_ice))/mu0(ig))
2704            else
2705             ! only dust contributes
2706            fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_sol*mu0(ig)*
2707     &                    exp(-(tau(ig,iaer_dust_doubleq))/mu0(ig))
2708            endif ! of if (water)
2709           endif ! of if(mu0(ig).le.0.01)
2710         enddo
2711
2712           ! Density-scaled opacities
2713              do ig=1,ngrid
2714                dsodust(ig,:) =
2715     &                dsodust(ig,:)*tauscaling(ig)
2716                dsords(ig,:) =
2717     &                dsords(ig,:)*tauscaling(ig)
2718                dsotop(ig,:) =
2719     &                dsotop(ig,:)*tauscaling(ig)
2720              enddo
2721
2722           if(doubleq) then
2723              do ig=1,ngrid
2724         IF (sedimentation) THEN
2725                dqdustsurf(ig) =
2726     &                zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
2727                dndustsurf(ig) =
2728     &                zdqssed(ig,igcm_dust_number)*tauscaling(ig)
2729         ENDIF
2730                ndust(ig,:) =
2731     &                zq(ig,:,igcm_dust_number)*tauscaling(ig)
2732                qdust(ig,:) =
2733     &                zq(ig,:,igcm_dust_mass)*tauscaling(ig)
2734              enddo
2735              if (scavenging) then
2736                do ig=1,ngrid
2737         IF (sedimentation) THEN
2738                  dqdustsurf(ig) = dqdustsurf(ig) +
2739     &                     zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
2740                  dndustsurf(ig) = dndustsurf(ig) +
2741     &                     zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
2742         ENDIF
2743                  nccn(ig,:) =
2744     &                     zq(ig,:,igcm_ccn_number)*tauscaling(ig)
2745                  qccn(ig,:) =
2746     &                     zq(ig,:,igcm_ccn_mass)*tauscaling(ig)
2747                enddo
2748              endif
2749           endif ! of (doubleq)
2750
2751           if (rdstorm) then   ! diagnostics of stormdust tendancies for 1D and 3D
2752              mstormdtot(:)=0
2753              mdusttot(:)=0
2754              qdusttotal(:,:)=0
2755              do ig=1,ngrid
2756                rdsdqdustsurf(ig) =
2757     &                zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig)
2758                rdsdndustsurf(ig) =
2759     &                zdqssed(ig,igcm_stormdust_number)*tauscaling(ig)
2760                rdsndust(ig,:) =
2761     &                pq(ig,:,igcm_stormdust_number)*tauscaling(ig)
2762                rdsqdust(ig,:) =
2763     &                pq(ig,:,igcm_stormdust_mass)*tauscaling(ig)
2764                do l=1,nlayer
2765                    mstormdtot(ig) = mstormdtot(ig) +
2766     &                      zq(ig,l,igcm_stormdust_mass) *
2767     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
2768                    mdusttot(ig) = mdusttot(ig) +
2769     &                      zq(ig,l,igcm_dust_mass) *
2770     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
2771                    qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust
2772                enddo
2773              enddo
2774           endif !(rdstorm)
2775
2776           if (water) then
2777             mtot(:)=0
2778             icetot(:)=0
2779             rave(:)=0
2780             tauTES(:)=0
2781
2782             IF (hdo) then
2783                 mtotD(:)=0
2784                 icetotD(:)=0
2785             ENDIF !hdo
2786
2787             do ig=1,ngrid
2788               do l=1,nlayer
2789                 mtot(ig) = mtot(ig) +
2790     &                      zq(ig,l,igcm_h2o_vap) *
2791     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
2792                 icetot(ig) = icetot(ig) +
2793     &                        zq(ig,l,igcm_h2o_ice) *
2794     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
2795                 IF (hdo) then
2796                 mtotD(ig) = mtotD(ig) +
2797     &                      zq(ig,l,igcm_hdo_vap) *
2798     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
2799                 icetotD(ig) = icetotD(ig) +
2800     &                        zq(ig,l,igcm_hdo_ice) *
2801     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
2802                 ENDIF !hdo
2803
2804c                Computing abs optical depth at 825 cm-1 in each
2805c                  layer to simulate NEW TES retrieval
2806                 Qabsice = min(
2807     &             max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2
2808     &                        )
2809                 opTES(ig,l)= 0.75 * Qabsice *
2810     &             zq(ig,l,igcm_h2o_ice) *
2811     &             (zplev(ig,l) - zplev(ig,l+1)) / g
2812     &             / (rho_ice * rice(ig,l) * (1.+nuice_ref))
2813                 tauTES(ig)=tauTES(ig)+ opTES(ig,l)
2814               enddo
2815c              rave(ig)=rave(ig)/max(icetot(ig),1.e-30)       ! mass weight
2816c               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
2817             enddo
2818             call watersat(ngrid*nlayer,zt,zplay,zqsat)
2819             satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
2820
2821             if (scavenging) then
2822               Nccntot(:)= 0
2823               Mccntot(:)= 0
2824               rave(:)=0
2825               do ig=1,ngrid
2826                 do l=1,nlayer
2827                    Nccntot(ig) = Nccntot(ig) +
2828     &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
2829     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
2830                    Mccntot(ig) = Mccntot(ig) +
2831     &              zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
2832     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
2833cccc Column integrated effective ice radius
2834cccc is weighted by total ice surface area (BETTER than total ice mass)
2835                    rave(ig) = rave(ig) +
2836     &                      tauscaling(ig) *
2837     &                      zq(ig,l,igcm_ccn_number) *
2838     &                      (zplev(ig,l) - zplev(ig,l+1)) / g *
2839     &                      rice(ig,l) * rice(ig,l)*  (1.+nuice_ref)
2840                 enddo
2841               rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
2842     &                  /max(pi*rave(ig),1.e-30) ! surface weight
2843               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
2844               enddo
2845             else ! of if (scavenging)
2846               rave(:)=0
2847               do ig=1,ngrid
2848                 do l=1,nlayer
2849                 rave(ig) = rave(ig) +
2850     &                      zq(ig,l,igcm_h2o_ice) *
2851     &                      (zplev(ig,l) - zplev(ig,l+1)) / g *
2852     &                      rice(ig,l) * (1.+nuice_ref)
2853                 enddo
2854                 rave(ig) = max(rave(ig) /
2855     &             max(icetot(ig),1.e-30),1.e-30) ! mass weight
2856               enddo
2857             endif ! of if (scavenging)
2858
2859           !Alternative A. Pottier weighting
2860           rave2(:) = 0.
2861           totrave2(:) = 0.
2862           do ig=1,ngrid
2863              do l=1,nlayer
2864              rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l)
2865              totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice)
2866              end do
2867              rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30)
2868           end do
2869
2870           endif ! of if (water)
2871
2872          if (co2clouds) then
2873            mtotco2(1:ngrid) = 0.
2874            icetotco2(1:ngrid) = 0.
2875            vaptotco2(1:ngrid) = 0.
2876            do ig=1,ngrid
2877              do l=1,nlayer
2878                vaptotco2(ig) = vaptotco2(ig) +
2879     &                          zq(ig,l,igcm_co2) *
2880     &                          (zplev(ig,l) - zplev(ig,l+1)) / g
2881                icetotco2(ig) = icetot(ig) +
2882     &                          zq(ig,l,igcm_co2_ice) *
2883     &                          (zplev(ig,l) - zplev(ig,l+1)) / g
2884              end do
2885              mtotco2(ig) = icetotco2(ig) + vaptotco2(ig)
2886            end do
2887          end if
2888
2889#ifndef MESOSCALE
2890c        -----------------------------------------------------------------
2891c        WSTATS: Saving statistics
2892c        -----------------------------------------------------------------
2893c        ("stats" stores and accumulates key variables in file "stats.nc"
2894c        which can later be used to make the statistic files of the run:
2895c        if flag "callstats" from callphys.def is .true.)
2896
2897        call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
2898        call wstats(ngrid,"tsurf","Surface temperature","K",2
2899     &      ,tsurf(:,iflat))
2900        call wstats(ngrid,"co2ice","CO2 ice cover",
2901     &                "kg.m-2",2,qsurf(:,igcm_co2,iflat))
2902        call wstats(ngrid,"watercap","H2O ice cover",
2903     &                "kg.m-2",2,watercap(:,iflat))
2904        call wstats(ngrid,"tau_pref_scenario",
2905     &              "prescribed visible dod at 610 Pa","NU",
2906     &                2,tau_pref_scenario)
2907        call wstats(ngrid,"tau_pref_gcm",
2908     &              "visible dod at 610 Pa in the GCM","NU",
2909     &                2,tau_pref_gcm)
2910        call wstats(ngrid,"fluxsurf_lw",
2911     &                "Thermal IR radiative flux to surface","W.m-2",2,
2912     &                fluxsurf_lw(:,iflat))
2913        call wstats(ngrid,"fluxsurf_dn_sw",
2914     &        "Incoming Solar radiative flux to surface","W.m-2",2,
2915     &                fluxsurf_dn_sw_tot(:,iflat))
2916        call wstats(ngrid,"fluxsurf_up_sw",
2917     &        "Reflected Solar radiative flux from surface","W.m-2",2,
2918     &                fluxsurf_up_sw_tot)
2919        call wstats(ngrid,"fluxtop_lw",
2920     &                "Thermal IR radiative flux to space","W.m-2",2,
2921     &                fluxtop_lw)
2922        call wstats(ngrid,"fluxtop_dn_sw",
2923     &        "Incoming Solar radiative flux from space","W.m-2",2,
2924     &                fluxtop_dn_sw_tot)
2925        call wstats(ngrid,"fluxtop_up_sw",
2926     &        "Outgoing Solar radiative flux to space","W.m-2",2,
2927     &                fluxtop_up_sw_tot)
2928        call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
2929        call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
2930        call wstats(ngrid,"v","Meridional (North-South) wind",
2931     &                "m.s-1",3,zv)
2932        call wstats(ngrid,"w","Vertical (down-up) wind",
2933     &                "m.s-1",3,pw)
2934        call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
2935        call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
2936          call wstats(ngrid,"q2",
2937     &                "Boundary layer eddy kinetic energy",
2938     &                "m2.s-2",3,q2)
2939          call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
2940     &                emis(:,iflat))
2941          call wstats(ngrid,"fluxsurf_dir_dn_sw",
2942     &                "Direct incoming SW flux at surface",
2943     &                "W.m-2",2,fluxsurf_dir_dn_sw)
2944
2945          if (calltherm) then
2946            call wstats(ngrid,"zmax_th","Height of thermals",
2947     &                "m",2,zmax_th)
2948            call wstats(ngrid,"hfmax_th","Max thermals heat flux",
2949     &                "K.m/s",2,hfmax_th)
2950            call wstats(ngrid,"wstar",
2951     &                "Max vertical velocity in thermals",
2952     &                "m/s",2,wstar)
2953          endif
2954
2955             if (water) then
2956               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2957     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2958               call wstats(ngrid,"vmr_h2ovap",
2959     &                    "H2O vapor volume mixing ratio","mol/mol",
2960     &                    3,vmr)
2961               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2962     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2963               call wstats(ngrid,"vmr_h2oice",
2964     &                    "H2O ice volume mixing ratio","mol/mol",
2965     &                    3,vmr)
2966              ! also store vmr_ice*rice for better diagnostics of rice
2967               vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
2968     &                               rice(1:ngrid,1:nlayer)
2969               call wstats(ngrid,"vmr_h2oice_rice",
2970     &                "H2O ice mixing ratio times ice particule size",
2971     &                    "(mol/mol)*m",
2972     &                    3,vmr)
2973               vmr=zqsat(1:ngrid,1:nlayer)
2974     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2975               call wstats(ngrid,"vmr_h2osat",
2976     &                    "saturation volume mixing ratio","mol/mol",
2977     &                    3,vmr)
2978               call wstats(ngrid,"h2o_ice_s",
2979     &                    "surface h2o_ice","kg/m2",
2980     &                    2,qsurf(1,igcm_h2o_ice,iflat))
2981               call wstats(ngrid,'albedo',
2982     &                         'albedo',
2983     &                         '',2,albedo(1,1,iflat))
2984               call wstats(ngrid,"mtot",
2985     &                    "total mass of water vapor","kg/m2",
2986     &                    2,mtot)
2987               call wstats(ngrid,"icetot",
2988     &                    "total mass of water ice","kg/m2",
2989     &                    2,icetot)
2990               call wstats(ngrid,"reffice",
2991     &                    "Mean reff","m",
2992     &                    2,rave)
2993              call wstats(ngrid,"Nccntot",
2994     &                    "condensation nuclei","Nbr/m2",
2995     &                    2,Nccntot)
2996              call wstats(ngrid,"Mccntot",
2997     &                    "condensation nuclei mass","kg/m2",
2998     &                    2,Mccntot)
2999              call wstats(ngrid,"rice",
3000     &                    "Ice particle size","m",
3001     &                    3,rice)
3002               if (.not.activice) then
3003                 call wstats(ngrid,"tauTESap",
3004     &                      "tau abs 825 cm-1","",
3005     &                      2,tauTES)
3006               else
3007                 call wstats(ngrid,'tauTES',
3008     &                   'tau abs 825 cm-1',
3009     &                  '',2,taucloudtes)
3010               endif
3011
3012             endif ! of if (water)
3013
3014             if (co2clouds) then
3015               call wstats(ngrid,"mtotco2",
3016     &                    "total mass atm of co2","kg/m2",
3017     &                    2,mtotco2)
3018               call wstats(ngrid,"icetotco2",
3019     &                    "total mass atm of co2 ice","kg/m2",
3020     &                    2,icetotco2)
3021               call wstats(ngrid,"vaptotco2",
3022     &                    "total mass atm of co2 vapor","kg/m2",
3023     &                    2,vaptotco2)
3024             end if
3025
3026
3027           if (dustbin.ne.0) then
3028
3029             call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
3030
3031             if (doubleq) then
3032               call wstats(ngrid,'dqsdust',
3033     &                        'deposited surface dust mass',
3034     &                        'kg.m-2.s-1',2,dqdustsurf)
3035               call wstats(ngrid,'dqndust',
3036     &                        'deposited surface dust number',
3037     &                        'number.m-2.s-1',2,dndustsurf)
3038               call wstats(ngrid,'reffdust','reffdust',
3039     &                        'm',3,rdust*ref_r0)
3040               call wstats(ngrid,'dustq','Dust mass mr',
3041     &                        'kg/kg',3,qdust)
3042               call wstats(ngrid,'dustN','Dust number',
3043     &                        'part/kg',3,ndust)
3044               if (rdstorm) then
3045                 call wstats(ngrid,'reffstormdust','reffdust',
3046     &                          'm',3,rstormdust*ref_r0)
3047                 call wstats(ngrid,'rdsdustq','Dust mass mr',
3048     &                          'kg/kg',3,rdsqdust)
3049                 call wstats(ngrid,'rdsdustN','Dust number',
3050     &                        'part/kg',3,rdsndust)
3051               end if
3052             else
3053               do iq=1,dustbin
3054                 write(str2(1:2),'(i2.2)') iq
3055                 call wstats(ngrid,'q'//str2,'mix. ratio',
3056     &                         'kg/kg',3,zq(1,1,iq))
3057                 call wstats(ngrid,'qsurf'//str2,'qsurf',
3058     &                         'kg.m-2',2,qsurf(1,iq,iflat))
3059               end do
3060             endif ! (doubleq)
3061
3062             if (scavenging) then
3063               call wstats(ngrid,'ccnq','CCN mass mr',
3064     &                        'kg/kg',3,qccn)
3065               call wstats(ngrid,'ccnN','CCN number',
3066     &                        'part/kg',3,nccn)
3067             endif ! (scavenging)
3068
3069           endif ! (dustbin.ne.0)
3070
3071           if (photochem) then
3072              do iq=1,nq
3073                 if (noms(iq) .ne. "dust_mass" .and.
3074     $               noms(iq) .ne. "dust_number" .and.
3075     $               noms(iq) .ne. "ccn_mass" .and.
3076     $               noms(iq) .ne. "ccn_number" .and.
3077     $               noms(iq) .ne. "ccnco2_mass" .and.
3078     $               noms(iq) .ne. "ccnco2_number" .and.
3079     $               noms(iq) .ne. "stormdust_mass" .and.
3080     $               noms(iq) .ne. "stormdust_number" .and.
3081     $               noms(iq) .ne. "topdust_mass" .and.
3082     $               noms(iq) .ne. "topdust_number") then
3083!                   volume mixing ratio
3084
3085                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
3086     &                            *mmean(1:ngrid,1:nlayer)/mmol(iq)
3087
3088                    call wstats(ngrid,"vmr_"//trim(noms(iq)),
3089     $                        "Volume mixing ratio","mol/mol",3,vmr)
3090                    if ((noms(iq).eq."o")
3091     $             .or. (noms(iq).eq."co2")
3092     $             .or. (noms(iq).eq."o3")
3093     $             .or. (noms(iq).eq."ar")
3094     $             .or. (noms(iq).eq."o2")
3095     $             .or. (noms(iq).eq."h2o_vap") ) then
3096                      call write_output("vmr_"//trim(noms(iq)),
3097     $                         "Volume mixing ratio","mol/mol",vmr(:,:))
3098                    end if
3099
3100!                   number density (molecule.cm-3)
3101
3102                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
3103     &                          *rho(1:ngrid,1:nlayer)*n_avog/
3104     &                           (1000*mmol(iq))
3105
3106                   call wstats(ngrid,"num_"//trim(noms(iq)),
3107     $                   "Number density","cm-3",3,rhopart)
3108                   call write_output("num_"//trim(noms(iq)),
3109     $                  "Number density","cm-3",rhopart(:,:))
3110
3111!                   vertical column (molecule.cm-2)
3112
3113                    do ig = 1,ngrid
3114                       colden(ig,iq) = 0.
3115                    end do
3116                    do l=1,nlayer
3117                       do ig=1,ngrid
3118                          colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
3119     $                                   *(zplev(ig,l)-zplev(ig,l+1))
3120     $                                   *6.022e22/(mmol(iq)*g)
3121                       end do
3122                    end do
3123
3124                    call wstats(ngrid,"c_"//trim(noms(iq)),
3125     $                          "column","mol cm-2",2,colden(1,iq))
3126                    call write_output("c_"//trim(noms(iq)),
3127     $                          "column","mol cm-2",colden(:,iq))
3128
3129!                   global mass (g)
3130
3131                    call planetwide_sumval(colden(:,iq)/6.022e23
3132     $                            *mmol(iq)*1.e4*cell_area(:),mass(iq))
3133
3134                    call write_output("mass_"//trim(noms(iq)),
3135     $                              "global mass","g",mass(iq))
3136
3137                 end if ! of if (noms(iq) .ne. "dust_mass" ...)
3138              end do ! of do iq=1,nq
3139           end if ! of if (photochem)
3140
3141           IF(lastcall.and.callstats) THEN
3142             write (*,*) "Writing stats..."
3143             call mkstats(ierr)
3144           ENDIF
3145
3146c        (Store EOF for Mars Climate database software)
3147         IF (calleofdump) THEN
3148            CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
3149         ENDIF
3150#endif
3151!endif of ifndef MESOSCALE
3152
3153#ifdef MESOSCALE
3154
3155      !! see comm_wrf.
3156      !! not needed when an array is already in a shared module.
3157      !! --> example : hfmax_th, zmax_th
3158
3159      CALL allocate_comm_wrf(ngrid,nlayer)
3160
3161      !state  real  HR_SW     ikj   misc  1  -  h  "HR_SW"     "HEATING RATE SW"                 "K/s"
3162      comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
3163      !state  real  HR_LW     ikj   misc  1  -  h  "HR_LW"     "HEATING RATE LW"                 "K/s"
3164      comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
3165      !state  real  SWDOWNZ    ij   misc  1  -  h  "SWDOWNZ"   "DOWNWARD SW FLUX AT SURFACE"     "W m-2"
3166      comm_SWDOWNZ(1:ngrid) = fluxsurf_dn_sw_tot(1:ngrid)
3167      !state  real  TAU_DUST   ij   misc  1  -  h  "TAU_DUST"  "REFERENCE VISIBLE DUST OPACITY"  ""
3168      comm_TAU_DUST(1:ngrid) = tau_pref_gcm(1:ngrid)
3169      !state  real  RDUST     ikj   misc  1  -  h  "RDUST"     "DUST RADIUS"                     "m"
3170      comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
3171      !state  real  QSURFDUST  ij   misc  1  -  h  "QSURFDUST" "DUST MASS AT SURFACE"            "kg m-2"
3172      IF (igcm_dust_mass .ne. 0) THEN
3173        comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
3174      ELSE
3175        comm_QSURFDUST(1:ngrid) = 0.
3176      ENDIF
3177      !state  real  MTOT       ij   misc  1  -  h  "MTOT"      "TOTAL MASS WATER VAPOR in pmic"  "pmic"
3178      comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
3179      !state  real  ICETOT     ij   misc  1  -  h  "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"
3180      comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
3181      !state  real  VMR_ICE   ikj   misc  1  -  h  "VMR_ICE"   "VOL. MIXING RATIO ICE"           "ppm"
3182      IF (igcm_h2o_ice .ne. 0) THEN
3183        comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
3184     .    * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
3185     .    * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
3186      ELSE
3187        comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
3188      ENDIF
3189      !state  real  TAU_ICE    ij   misc  1  -  h  "TAU_ICE"   "CLOUD OD at 825 cm-1 TES"        ""
3190      if (activice) then
3191        comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
3192      else
3193        comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
3194      endif
3195      !state  real  RICE      ikj   misc  1  -  h  "RICE"      "ICE RADIUS"                      "m"
3196      comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
3197
3198      !! calculate sensible heat flux in W/m2 for outputs
3199      !! -- the one computed in vdifc is not the real one
3200      !! -- vdifc must have been called
3201      if (.not.callrichsl) then
3202        sensibFlux(1:ngrid) = zflubid(1:ngrid)
3203     .         - capcal(1:ngrid)*zdtsdif(1:ngrid)
3204      else
3205        sensibFlux(1:ngrid) =
3206     &   (pplay(1:ngrid,1)/(r*pt(1:ngrid,1)))*cpp
3207     &   *sqrt(pu(1:ngrid,1)*pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
3208     &         +(log(1.+0.7*wstar(1:ngrid) + 2.3*wstar(1:ngrid)**2))**2)
3209     &   *zcdh(1:ngrid)*(tsurf(1:ngrid)-zh(1:ngrid,1))
3210      endif
3211
3212#else
3213#ifndef MESOINI
3214
3215c        ==========================================================
3216c        WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
3217c          any variable for diagnostic (output with period
3218c          "ecritphy", set in "run.def")
3219c        ==========================================================
3220c        WRITEDIAGFI can ALSO be called from any other subroutines
3221c        for any variables !!
3222         call write_output("emis","Surface emissivity","",
3223     &                  emis(:,iflat))
3224         do islope=1,nslope
3225           write(str2(1:2),'(i2.2)') islope
3226           call write_output("emis_slope"//str2,
3227     &    "Surface emissivity","",emis(:,islope))
3228         ENDDO
3229         call write_output("zzlay","Midlayer altitude",
3230     &                    "m",zzlay(:,:))
3231         call write_output("zzlev","Interlayer altitude",
3232     &                    "m",zzlev(:,1:nlayer))
3233         call write_output("pphi","Geopotential","m2s-2",
3234     &                    pphi(:,:))
3235         call write_output("phisfi","Surface geopotential",
3236     &                    "m2s-2",phisfi(:))
3237         if (grid_type == regular_lonlat) then
3238           call write_output("area","Mesh area","m2",
3239     &                       cell_area_for_lonlat_outputs)
3240         else ! unstructured grid (e.g. dynamico)
3241           call write_output("area","Mesh area","m2",cell_area)
3242         endif
3243         call write_output("tsurf","Surface temperature","K",
3244     &                  tsurf(:,iflat))
3245         do islope=1,nslope
3246          write(str2(1:2),'(i2.2)') islope
3247         call write_output("tsurf_slope"//str2,
3248     &            "Surface temperature","K",
3249     &                  tsurf(:,islope))
3250         ENDDO
3251         call write_output("ps","surface pressure","Pa",ps(:))
3252         call write_output("co2ice","co2 ice thickness"
3253     &                              ,"kg.m-2",qsurf(:,igcm_co2,iflat))
3254         do islope=1,nslope
3255          write(str2(1:2),'(i2.2)') islope
3256         call write_output("co2ice_slope"//str2,"co2 ice thickness"
3257     &              ,"kg.m-2",qsurf(:,igcm_co2,islope))
3258         ENDDO
3259         call write_output("watercap","Perennial water ice thickness"
3260     &         ,"kg.m-2",watercap(:,iflat))
3261         do islope=1,nslope
3262          write(str2(1:2),'(i2.2)') islope
3263         call write_output("watercap_slope"//str2,
3264     &         "Perennial water ice thickness"
3265     &         ,"kg.m-2",watercap(:,islope))
3266         ENDDO
3267         call write_output("perennial_co2ice",
3268     &         "Perennial co2 ice thickness","kg.m-2",
3269     &         perennial_co2ice(:,iflat))
3270         do islope=1,nslope
3271          write(str2(1:2),'(i2.2)') islope
3272         call write_output("perennial_co2ice_slope"//str2,
3273     &         "Perennial co2 ice thickness"
3274     &         ,"kg.m-2",perennial_co2ice(:,islope))
3275         ENDDO
3276         call write_output("temp_layer1","temperature in layer 1",
3277     &                  "K",zt(:,1))
3278         call write_output("temp7","temperature in layer 7",
3279     &                  "K",zt(:,7))
3280         call write_output("fluxsurf_lw","fluxsurf_lw","W.m-2",
3281     &                  fluxsurf_lw(:,iflat))
3282         do islope=1,nslope
3283          write(str2(1:2),'(i2.2)') islope
3284         call write_output("fluxsurf_lw_slope"//str2,
3285     &              "fluxsurf_lw","W.m-2",
3286     &                  fluxsurf_lw(:,islope))
3287         ENDDO
3288         call write_output("fluxsurf_dn_sw","fluxsurf_dn_sw",
3289     &                  "W.m-2",fluxsurf_dn_sw_tot(:,iflat))
3290         do islope=1,nslope
3291          write(str2(1:2),'(i2.2)') islope
3292         call write_output("fluxsurf_dn_sw_slope"//str2,
3293     &                  "fluxsurf_dn_sw",
3294     &                  "W.m-2",fluxsurf_dn_sw_tot(:,islope))
3295         ENDDO
3296         call write_output("fluxtop_dn_sw","fluxtop_dn_sw",
3297     &                  "W.m-2",fluxtop_dn_sw(:,1) + fluxtop_dn_sw(:,2))
3298         call write_output("fluxtop_lw","fluxtop_lw","W.m-2",
3299     &                  fluxtop_lw(:))
3300         call write_output("fluxtop_up_sw","fluxtop_up_sw","W.m-2",
3301     &                  fluxtop_up_sw_tot(:))
3302         call write_output("temp","temperature","K",zt(:,:))
3303         call write_output("Sols","Time","sols",zday)
3304         call write_output("Ls","Solar longitude","deg",
3305     &                    zls*180./pi)
3306         call write_output("u","Zonal wind","m.s-1",zu(:,:))
3307         call write_output("v","Meridional wind","m.s-1",zv(:,:))
3308         call write_output("w","Vertical wind","m.s-1",pw(:,:))
3309         call write_output("rho","density","kg.m-3",rho(:,:))
3310         call write_output("pressure","Pressure","Pa",zplay(:,:))
3311         call write_output("zplev","Interlayer pressure","Pa",
3312     &                     zplev(:,1:nlayer))
3313        call write_output('sw_htrt','sw heat. rate',
3314     &                   'K/s',zdtsw(:,:))
3315        call write_output('lw_htrt','lw heat. rate',
3316     &                   'K/s',zdtlw(:,:))
3317        call write_output("local_time","Local time",
3318     &                   'sol',local_time(:))
3319        if (water) then
3320            if (.not.activice) then
3321               CALL write_output('tauTESap',
3322     &                         'tau abs 825 cm-1',
3323     &                         '',tauTES(:))
3324             else
3325               CALL write_output('tauTES',
3326     &                         'tau abs 825 cm-1',
3327     &                         '',taucloudtes(:))
3328             endif
3329        endif ! of if (water)
3330#else
3331      !!! this is to ensure correct initialisation of mesoscale model
3332      call write_output("tsurf","Surface temperature","K",
3333     &                 tsurf(:,iflat))
3334      call write_output("ps","surface pressure","Pa",ps(:))
3335      call write_output("co2ice","co2 ice thickness","kg.m-2",
3336     &                 qsurf(:,igcm_co2,iflat))
3337      call write_output("temp","temperature","K",zt(:,:))
3338      call write_output("u","Zonal wind","m.s-1",zu(:,:))
3339      call write_output("v","Meridional wind","m.s-1",zv(:,:))
3340      call write_output("emis","Surface emissivity","",
3341     &                 emis(:,iflat))
3342      call write_output("tsoil","Soil temperature",
3343     &                 "K",tsoil(:,:,iflat))
3344      call write_output("inertiedat","Soil inertia",
3345     &                 "K",inertiedat(:,:))
3346#endif
3347
3348c        ----------------------------------------------------------
3349c        Outputs of the CO2 cycle
3350c        ----------------------------------------------------------
3351      if (igcm_co2.ne.0) then
3352        call write_output("co2","co2 mass mixing ratio",
3353     &                   "kg.kg-1",zq(:,:,igcm_co2))
3354
3355        if (co2clouds) then
3356          call write_output('ccnqco2','CCNco2 mmr',
3357     &                     'kg.kg-1',zq(:,:,igcm_ccnco2_mass))
3358
3359          call write_output('ccnNco2','CCNco2 number',
3360     &                     'part.kg-1',zq(:,:,igcm_ccnco2_number))
3361
3362          call write_output('co2_ice','co2_ice mmr in atm',
3363     &                     'kg.kg-1',zq(:,:,igcm_co2_ice))
3364
3365         call write_output("mtotco2","total mass atm of co2",
3366     &                    "kg.m-2",mtotco2(:))
3367         call write_output("icetotco2","total mass atm of co2 ice",
3368     &                    "kg.m-2", icetotco2(:))
3369         call write_output("vaptotco2","total mass atm of co2 "//
3370     &                    "vapor","kg.m-2", vaptotco2(:))
3371         if (co2useh2o) then
3372           call write_output('ccnqco2_h2o_m_ice',
3373     &                      'CCNco2_h2o_mass_ice mmr',
3374     &                    'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ice))
3375
3376           call write_output('ccnqco2_h2o_m_ccn',
3377     &                      'CCNco2_h2o_mass_ccn mmr',
3378     &                   'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ccn))
3379
3380           call write_output('ccnNco2_h2o','CCNco2_h2o number',
3381     &                   'part.kg-1',zq(:,:,igcm_ccnco2_h2o_number))
3382         end if
3383
3384           if (meteo_flux) then
3385          call write_output('ccnqco2_meteor','CCNco2_meteor mmr',
3386     &                  'kg.kg-1',zq(:,:,igcm_ccnco2_meteor_mass))
3387
3388        call write_output('ccnNco2_meteor','CCNco2_meteor number',
3389     &                'part.kg-1',zq(:,:,igcm_ccnco2_meteor_number))
3390           end if
3391
3392        end if ! of if (co2clouds)
3393      end if ! of if (igcm_co2.ne.0)
3394
3395      ! Output He tracer, if there is one
3396      if (igcm_he.ne.0) then
3397        call write_output("he","helium mass mixing ratio",
3398     &                   "kg/kg",zq(:,:,igcm_he))
3399        vmr = zq(1:ngrid,1:nlayer,igcm_he)
3400     &        * mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
3401        call write_output('vmr_he','helium vol. mixing ratio',
3402     &                   'mol/mol',vmr(:,:))
3403      end if
3404
3405c        ----------------------------------------------------------
3406c        Outputs of the water cycle
3407c        ----------------------------------------------------------
3408        if (water) then
3409#ifdef MESOINI
3410            !!!! waterice = q01, voir readmeteo.F90
3411          call write_output('q01',noms(igcm_h2o_ice),
3412     &                     'kg/kg',
3413     &                     zq(:,:,igcm_h2o_ice))
3414            !!!! watervapor = q02, voir readmeteo.F90
3415          call write_output('q02',noms(igcm_h2o_vap),
3416     &                     'kg/kg',
3417     &                     zq(:,:,igcm_h2o_vap))
3418            !!!! surface waterice qsurf02 (voir readmeteo)
3419          call write_output('qsurf02','surface tracer',
3420     &                     'kg.m-2',
3421     &                     qsurf(:,igcm_h2o_ice,iflat))
3422         do islope=1,nslope
3423          write(str2(1:2),'(i2.2)') islope
3424          call write_output('qsurf02_slope'//str2,
3425     &         'surface tracer','kg.m-2',
3426     &                     qsurf(:,igcm_h2o_ice,islope))
3427         ENDDO
3428#endif
3429          call write_output('mtot',
3430     &                     'total mass of water vapor',
3431     &                     'kg/m2',mtot(:))
3432          call write_output('icetot',
3433     &                     'total mass of water ice',
3434     &                     'kg/m2',icetot(:))
3435          vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice)
3436     &          * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
3437          call write_output('vmr_h2oice','h2o ice vmr',
3438     &                     'mol/mol',vmr(:,:))
3439          vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap)
3440     &          * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
3441          call write_output('vmr_h2ovap','h2o vap vmr',
3442     &                     'mol/mol',vmr(:,:))
3443          call write_output('reffice',
3444     &                     'Mean reff',
3445     &                     'm',rave(:))
3446          call write_output('h2o_ice','h2o_ice','kg/kg',
3447     &                     zq(:,:,igcm_h2o_ice))
3448          call write_output('h2o_vap','h2o_vap','kg/kg',
3449     &                     zq(:,:,igcm_h2o_vap))
3450
3451            if (hdo) then
3452            vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice)
3453     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice)
3454            CALL write_output('vmr_hdoice','hdo ice vmr',
3455     &                       'mol/mol',vmr(:,:))
3456            vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap)
3457     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap)
3458            CALL write_output('vmr_hdovap','hdo vap vmr',
3459     &                       'mol/mol',vmr(:,:))
3460            call write_output('hdo_ice','hdo_ice','kg/kg',
3461     &             zq(:,:,igcm_hdo_ice))
3462            call write_output('hdo_vap','hdo_vap','kg/kg',
3463     &             zq(:,:,igcm_hdo_vap))
3464
3465            CALL write_output('mtotD',
3466     &                       'total mass of HDO vapor',
3467     &                       'kg/m2',mtotD(:))
3468            CALL write_output('icetotD',
3469     &                       'total mass of HDO ice',
3470     &                       'kg/m2',icetotD(:))
3471
3472C           Calculation of the D/H ratio
3473            do l=1,nlayer
3474                do ig=1,ngrid
3475                if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
3476                    DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
3477     &              zq(ig,l,igcm_h2o_vap) )*1./(2.*155.76e-6)
3478                else
3479                    DoH_vap(ig,l) = 0.
3480                endif
3481                enddo
3482            enddo
3483
3484            do l=1,nlayer
3485                do ig=1,ngrid
3486                if (zq(ig,l,igcm_h2o_ice).gt.qperemin) then
3487                    DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/
3488     &                  zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6)
3489                else
3490                    DoH_ice(ig,l) = 0.
3491                endif
3492                enddo
3493            enddo
3494
3495            CALL write_output('DoH_vap',
3496     &                       'D/H ratio in vapor',
3497     &                       ' ',DoH_vap(:,:))
3498            CALL write_output('DoH_ice',
3499     &                       'D/H ratio in ice',
3500     &                       '',DoH_ice(:,:))
3501
3502            endif !hdo
3503
3504!A. Pottier
3505!             CALL write_output('rmoym',
3506!     &                      'alternative reffice',
3507!     &                      'm',rave2(:))
3508            call write_output('h2o_saturation',
3509     &           'h2o vap saturation ratio','',satu(:,:))
3510            if (scavenging) then
3511              CALL write_output("Nccntot",
3512     &                    "condensation nuclei","Nbr/m2",
3513     &                    Nccntot(:))
3514              CALL write_output("Mccntot",
3515     &                    "mass condensation nuclei","kg/m2",
3516     &                    Mccntot(:))
3517            endif
3518            call write_output('rice','Water ice particle size',
3519     &                       'm',rice(:,:))
3520            call write_output('h2o_ice_s',
3521     &                       'surface h2o_ice',
3522     &            'kg.m-2',qsurf(:,igcm_h2o_ice,iflat))
3523         do islope=1,nslope
3524          write(str2(1:2),'(i2.2)') islope
3525            call write_output('h2o_ice_s_slope'//str2,
3526     &                       'surface h2o_ice',
3527     &             'kg.m-2',qsurf(:,igcm_h2o_ice,islope))
3528         ENDDO
3529            if (hdo) then
3530            call write_output('hdo_ice_s',
3531     &                       'surface hdo_ice',
3532     &           'kg.m-2',qsurf(:,igcm_hdo_ice,iflat))
3533         do islope=1,nslope
3534          write(str2(1:2),'(i2.2)') islope
3535            call write_output('hdo_ice_s_slope'//str2,
3536     &                       'surface hdo_ice',
3537     &           'kg.m-2',qsurf(:,igcm_hdo_ice,islope))
3538         ENDDO
3539
3540                do ig=1,ngrid
3541                if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then
3542           DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/
3543     &          qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6
3544                else
3545                    DoH_surf(ig) = 0.
3546                endif
3547                enddo
3548
3549            call write_output('DoH_surf',
3550     &                       'surface D/H',
3551     &                       '',DoH_surf(:))
3552            endif ! hdo
3553
3554            CALL write_output('albedo',
3555     &                         'albedo',
3556     &                         '',albedo(:,1,iflat))
3557         do islope=1,nslope
3558          write(str2(1:2),'(i2.2)') islope
3559            CALL write_output('albedo_slope'//str2,
3560     &                         'albedo',
3561     &                         '',albedo(:,1,islope))
3562         ENDDO
3563              if (surfaceice_tifeedback.or.poreice_tifeedback) then
3564                 call write_output("soiltemp",
3565     &                              "Soil temperature","K",
3566     &                              tsoil(:,:,iflat))
3567                 call write_output('soilti',
3568     &                       'Soil Thermal Inertia',
3569     &         'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat))
3570
3571         do islope=1,nslope
3572          write(str2(1:2),'(i2.2)') islope
3573                 call write_output('soilti_slope'//str2,
3574     &                       'Soil Thermal Inertia',
3575     &          'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,islope))
3576         ENDDO
3577              endif
3578!A. Pottier
3579          if (CLFvarying) then !AP14 nebulosity
3580            call write_output('totcloudfrac',
3581     &                       'Total cloud fraction',
3582     &                       ' ',totcloudfrac(:))
3583          end if !clf varying
3584        end if !(water)
3585
3586c        ----------------------------------------------------------
3587c        Outputs of the dust cycle
3588c        ----------------------------------------------------------
3589
3590      call write_output('tau_pref_scenario',
3591     &                 'Prescribed visible dust optical depth at 610Pa',
3592     &                 'NU',tau_pref_scenario(:))
3593
3594      call write_output('tau_pref_gcm',
3595     &                 'Visible dust optical depth at 610Pa in the GCM',
3596     &                 'NU',tau_pref_gcm(:))
3597
3598      if (reff_driven_IRtoVIS_scenario) then
3599             call write_output('IRtoVIScoef',
3600     &  'Conversion coeff for dust tau from abs9.3um to ext0.67um',
3601     &                        '/',IRtoVIScoef(:))
3602      endif
3603
3604      if (dustbin.ne.0) then
3605
3606#ifndef MESOINI
3607           if (doubleq) then
3608             call write_output('dqsdust',
3609     &                        'deposited surface dust mass',
3610     &                        'kg.m-2.s-1',dqdustsurf(:))
3611             call write_output('dqndust',
3612     &                        'deposited surface dust number',
3613     &                        'number.m-2.s-1',dndustsurf(:))
3614             call write_output('reffdust','reffdust',
3615     &                        'm',rdust(:,:)*ref_r0)
3616             call write_output('dustq','Dust mass mr',
3617     &                        'kg/kg',qdust(:,:))
3618             call write_output('dustN','Dust number',
3619     &                        'part/kg',ndust(:,:))
3620
3621         select case (trim(dustiropacity))
3622              case ("tes")
3623               call write_output('dsodust_TES',
3624     &  'density scaled extinction opacity of std dust at 9.3um(TES)',
3625     &         'm2.kg-1',dsodust(:,:))
3626               call write_output('dso_TES',
3627     &  'density scaled extinction opacity of all dust at 9.3um(TES)',
3628     &         'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
3629              case ("mcs")
3630               call write_output('dsodust',
3631     &  'density scaled extinction opacity of std dust at 21.6um(MCS)',
3632     &         'm2.kg-1',dsodust(:,:))
3633               call write_output('dso',
3634     &  'density scaled extinction opacity of all dust at 21.6um(MCS)',
3635     &         'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:))
3636             end select
3637           else ! (doubleq=.false.)
3638             do iq=1,dustbin
3639               write(str2(1:2),'(i2.2)') iq
3640               call write_output('q'//str2,'mix. ratio',
3641     &                         'kg/kg',zq(:,:,iq))
3642               call write_output('qsurf'//str2,'qsurf',
3643     &                         'kg.m-2',qsurf(:,iq,iflat))
3644         do islope=1,nslope
3645          write(str2(1:2),'(i2.2)') islope
3646               call write_output('qsurf_slope'//str2,'qsurf',
3647     &                         'kg.m-2',qsurf(:,iq,islope))
3648         ENDDO
3649             end do
3650           endif ! (doubleq)
3651
3652           if (rdstorm) then  ! writediagfi tendencies stormdust tracers
3653             call write_output('reffstormdust','reffstormdust',
3654     &                        'm',rstormdust(:,:)*ref_r0)
3655             call write_output('mstormdtot',
3656     &                        'total mass of stormdust only',
3657     &                        'kg.m-2',mstormdtot(:))
3658             call write_output('mdusttot',
3659     &                        'total mass of dust only',
3660     &                        'kg.m-2',mdusttot(:))
3661             call write_output('rdsdqsdust',
3662     &                        'deposited surface stormdust mass',
3663     &                        'kg.m-2.s-1',rdsdqdustsurf(:))
3664             call write_output('rdsdustq','storm Dust mass mr',
3665     &                        'kg/kg',rdsqdust(:,:))
3666             call write_output('rdsdustqmodel','storm Dust massmr',
3667     &                        'kg/kg',pq(:,:,igcm_stormdust_mass))
3668             call write_output('rdsdustN','storm Dust number',
3669     &                        'part/kg',rdsndust(:,:))
3670             call write_output("stormfract",
3671     &                "fraction of the mesh, with stormdust","none",
3672     &                                         totstormfract(:))
3673!             call write_output('qsurf',
3674!     &                 'stormdust injection',
3675!     &                 'kg.m-2',qsurf(:,igcm_stormdust_mass,iflat))
3676!         do islope=1,nslope
3677!          write(str2(1:2),'(i2.2)') islope
3678!             call write_output('qsurf_slope'//str2,
3679!     &                 'stormdust injection',
3680!     &          'kg.m-2',qsurf(:,igcm_stormdust_mass,islope))
3681!         ENDDO
3682!             call write_output('pdqsurf',
3683!     &                  'tendancy stormdust mass at surface',
3684!     &          'kg.m-2',dqsurf(:,igcm_stormdust_mass,iflat))
3685!         do islope=1,nslope
3686!          write(str2(1:2),'(i2.2)') islope
3687!             call write_output('pdqsurf_slope'//str2,
3688!     &                  'tendancy stormdust mass at surface',
3689!     &            'kg.m-2',dqsurf(:,igcm_stormdust_mass,islope))
3690!         ENDDO
3691             call write_output('wspeed_stormdust',
3692     &                         'vertical velocity of stormdust',
3693     &                         'm/s',wspeed(:,:))
3694             call write_output('zdqsed_dust_mass'
3695     &          ,'sedimentation tendency of background dust mmr'
3696     &          ,'kg/kg.s-1',
3697     &          zdqsed(:,:,igcm_dust_mass))
3698             call write_output('zdqssed_dust_mass'
3699     &          ,'sedimentation tendency of background dust on surface'
3700     &          ,'kg.m-2.s-1',
3701     &          zdqssed(:,igcm_dust_mass))
3702             call write_output('zdqsed_stormdust_mass'
3703     &          ,'sedimentation tendency of stormdust dust mmr'
3704     &          ,'kg/kg.s-1',
3705     &          zdqsed(:,:,igcm_stormdust_mass))
3706             call write_output('zdqsed_dust_number'
3707     &          ,'sedimentation tendency of background dust number'
3708     &          ,'nbr/kg.s-1',
3709     &          zdqsed(:,:,igcm_dust_number))
3710             call write_output('rdust','rdust',
3711     &                        'm',rdust(:,:))
3712             call write_output('rstormdust','rstormdust',
3713     &                        'm',rstormdust(:,:))
3714
3715         select case (trim(dustiropacity))
3716              case ("tes")
3717               call write_output('dsords_TES',
3718     &  'density scaled extinction opacity of stormdust at 9.3um(TES)',
3719     &         'm2.kg-1',dsords(:,:))
3720              case ("mcs")
3721               call write_output('dsords',
3722     &  'density scaled extinction opacity of stormdust at 21.6um(MCS)',
3723     &         'm2.kg-1',dsords(:,:))
3724             end select
3725           endif ! (rdstorm)
3726
3727           if (topflows) then
3728             call write_output('refftopdust',
3729     &                         'Topdust dust effective radius',
3730     &                         'm',rtopdust(:,:)*ref_r0)
3731             call write_output('topdustq','top Dust mass mr',
3732     &                        'kg/kg',pq(:,:,igcm_topdust_mass))
3733             call write_output('topdustN','top Dust number',
3734     &                        'part/kg',pq(:,:,igcm_topdust_number))
3735         select case (trim(dustiropacity))
3736              case ("tes")
3737               call write_output('dsotop_TES',
3738     &  'density scaled extinction opacity of topdust at 9.3um(TES)',
3739     &         'm2.kg-1',dsotop(:,:))
3740              case ("mcs")
3741               call write_output('dsotop',
3742     &  'density scaled extinction opacity of topdust at 21.6um(MCS)',
3743     &         'm2.kg-1',dsotop(:,:))
3744             end select
3745           endif ! (topflows)
3746
3747           if (dustscaling_mode==2) then
3748             call write_output("dust_rad_adjust",
3749     &            "radiative adjustment coefficient for dust",
3750     &                        "",dust_rad_adjust(:))
3751           endif
3752
3753!           if (scavenging) then ! these outputs should be in the scavenging routine
3754!             call write_output('ccnq','CCN mass mr',
3755!     &                        'kg/kg',qccn(:,:))
3756!             call write_output('ccnN','CCN number',
3757!     &                        'part/kg',nccn(:,:))
3758!             call write_output('surfccnq','Surf nuclei mass mr',
3759!     &                        'kg.m-2',qsurf(:,igcm_ccn_mass,iflat))
3760!         do islope=1,nslope
3761!          write(str2(1:2),'(i2.2)') islope
3762!             call write_output('surfccnq_slope'//str2,
3763!     &               'Surf nuclei mass mr',
3764!     &               'kg.m-2',qsurf(:,igcm_ccn_mass,islope))
3765!         ENDDO
3766!             call write_output('surfccnN','Surf nuclei number',
3767!     &                        'kg.m-2',qsurf(:,igcm_ccn_number,iflat))
3768!         do islope=1,nslope
3769!          write(str2(1:2),'(i2.2)') islope
3770!             call write_output('surfccnN_slope'//str2,
3771!     &                 'Surf nuclei number',
3772!     &                 'kg.m-2',qsurf(:,igcm_ccn_number,islope))
3773!         ENDDO
3774!           endif ! (scavenging)
3775
3776#else
3777!     !!! to initialize mesoscale we need scaled variables
3778!     !!! because this must correspond to starting point for tracers
3779!             call write_output('dustq','Dust mass mr',
3780!     &           'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
3781!             call write_output('dustN','Dust number',
3782!     &           'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
3783!             call write_output('ccn','Nuclei mass mr',
3784!     &           'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
3785!             call write_output('ccnN','Nuclei number',
3786!     &           'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
3787      if (freedust) then
3788             call write_output('dustq','Dust mass mr',
3789     &                        'kg/kg',qdust)
3790             call write_output('dustN','Dust number',
3791     &                        'part/kg',ndust)
3792             call write_output('ccn','CCN mass mr',
3793     &                        'kg/kg',qccn)
3794             call write_output('ccnN','CCN number',
3795     &                        'part/kg',nccn)
3796      else
3797             call write_output('dustq','Dust mass mr',
3798     &                        'kg/kg',pq(:,:,igcm_dust_mass))
3799             call write_output('dustN','Dust number',
3800     &                        'part/kg',pq(:,:,igcm_dust_number))
3801             call write_output('ccn','Nuclei mass mr',
3802     &                        'kg/kg',pq(:,:,igcm_ccn_mass))
3803             call write_output('ccnN','Nuclei number',
3804     &                        'part/kg',pq(:,:,igcm_ccn_number))
3805      endif
3806#endif
3807
3808         end if  ! (dustbin.ne.0)
3809
3810c        ----------------------------------------------------------
3811c        Thermospheric outputs
3812c        ----------------------------------------------------------
3813         if(callthermos) then
3814
3815            call write_output("quv","UV heating","K/s",
3816     $           zdteuv(:,:))
3817            call write_output("cond","Thermal conduction","K/s",
3818     $           zdtconduc(:,:))
3819
3820            !H, H2 and D escape fluxes
3821
3822            call write_output("PhiH","H escape flux","s-1",
3823     $           PhiEscH)
3824            call write_output("PhiH2","H2 escape flux","s-1",
3825     $           PhiEscH2)
3826            call write_output("PhiD","D escape flux","s-1",
3827     $           PhiEscD)
3828
3829         endif  !(callthermos)
3830
3831            call write_output("q15um","15 um cooling","K/s",
3832     $           zdtnlte(:,:))
3833            call write_output("qnir","NIR heating","K/s",
3834     $           zdtnirco2(:,:))
3835
3836c        ----------------------------------------------------------
3837c        ----------------------------------------------------------
3838c        PBL OUTPUS
3839c        ----------------------------------------------------------
3840c        ----------------------------------------------------------
3841
3842c        ----------------------------------------------------------
3843c        Outputs of thermals
3844c        ----------------------------------------------------------
3845      if (calltherm) then
3846        call write_output('lmax_th',
3847     &              'index of vertical extension of thermals',
3848     &                   'grid level',lmax_th_out(:))
3849        call write_output('zmax_th',
3850     &                   'vertical extension of thermals','m',
3851     &                    zmax_th(:))
3852        call write_output('hfmax_th',
3853     &                   'maximum heat flux in thermals','K.m/s',
3854     &                   hfmax_th(:))
3855        call write_output('wstar',
3856     &                   'maximum thermals vertical velocity','m/s',
3857     &                   wstar(:))
3858      end if
3859
3860c        ----------------------------------------------------------
3861c        ----------------------------------------------------------
3862c        END OF PBL OUTPUS
3863c        ----------------------------------------------------------
3864c        ----------------------------------------------------------
3865
3866c        ----------------------------------------------------------
3867c        Output in netcdf file "diagsoil.nc" for subterranean
3868c          variables (output every "ecritphy", as for writediagfi)
3869c        ----------------------------------------------------------
3870         ! Write soil temperature
3871        call write_output("soiltemp","Soil temperature","K",
3872     &                     tsoil(:,:,iflat))
3873         do islope=1,nslope
3874          write(str2(1:2),'(i2.2)') islope
3875        call write_output("soiltemp_slope"//str2,
3876     &                     "Soil temperature","K",
3877     &                     tsoil(:,:,islope))
3878         ENDDO
3879
3880!PREVIOUSLY IN 1D ONLY
3881         call write_output("dtrad","rad. heat. rate",
3882     &              "K.s-1",dtrad(:,:))
3883
3884           if (rdstorm) then
3885             call write_output('aerosol_dust','opacity of env. dust',''
3886     &                        ,aerosol(:,:,iaer_dust_doubleq))
3887             call write_output('aerosol_stormdust',
3888     &                         'opacity of storm dust',''
3889     &                        ,aerosol(:,:,iaer_stormdust_doubleq))
3890             call write_output('dqsdifdustq',
3891     &'tendency due to vertical diffusion of background dust on surface'
3892     &          ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,iflat))
3893         do islope=1,nslope
3894          write(str2(1:2),'(i2.2)') islope
3895             call write_output('dqsdifdustq_slope'//str2,
3896     &'tendency due to vertical diffusion of background dust on surface'
3897     &          ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,islope))
3898         ENDDO
3899             call write_output('dqsdifrdsq',
3900     & 'tendency due to vertical diffusion of stormdust on surface',
3901     &          'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,iflat))
3902         do islope=1,nslope
3903          write(str2(1:2),'(i2.2)') islope
3904             call write_output('dqsdifrdsq_slope'//str2,
3905     & 'tendency due to vertical diffusion of stormdust on surface',
3906     &          'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,islope))
3907         ENDDO
3908          endif !(rdstorm)
3909
3910         if(water) then
3911           if (.not.scavenging) then
3912        call write_output('zdqcloud_ice','cloud ice',
3913     &            'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice))
3914        call write_output('zdqcloud_vap','cloud vap',
3915     &            'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_vap))
3916        call write_output('zdqcloud','cloud',
3917     &            'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice)
3918     &                          +zdqcloud(:,:,igcm_h2o_vap))
3919        IF (hdo) THEN
3920        call write_output('zdqcloud_iceD','cloud ice hdo',
3921     &            'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_ice))
3922        call write_output('zdqcloud_vapD','cloud vap hdo',
3923     &            'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_vap))
3924        ENDIF ! hdo
3925           endif !not.scavenging
3926
3927! Output needed by the PEM
3928          DO ig = 1,ngrid
3929            ztmp1 =(1/m_co2 - 1/m_noco2)
3930            ztmp2=1/m_noco2
3931            pvap_surf(ig) = 1/(ztmp1*zq(ig,1,igcm_co2)+ztmp2)
3932     &      * zq(ig,1,igcm_h2o_vap)/(mmol(igcm_h2o_vap)*1.e-3)*ps(ig)
3933
3934            DO islope = 1,nslope
3935        ! Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
3936             rhowater_surf_sat(ig,islope)  =
3937     &         exp(beta_clap_h2o/tsurf(ig,islope)+alpha_clap_h2o)
3938     &         / tsurf(ig,islope)
3939     &         * mmol(igcm_h2o_vap)/(mugaz*r)
3940
3941             if(qsurf(ig,igcm_h2o_ice,islope).gt.(1.e-4)) then
3942           ! we consider to be at saturation above 1.e-4 kg.m-2
3943               rhowater_surf(ig,islope) = rhowater_surf_sat(ig,islope)
3944             else
3945          ! otherwise, use vapor partial pressure
3946              rhowater_surf(ig,islope) = pvap_surf(ig)
3947     &         / tsurf(ig,islope)
3948     &         * mmol(igcm_h2o_vap)/(mugaz*r)
3949             endif
3950             DO isoil = 1,nsoilmx
3951             rhowater_soil(ig,isoil,islope) =
3952     &         exp(beta_clap_h2o/tsoil(ig,isoil,islope)+alpha_clap_h2o)
3953     &         / tsoil(ig,isoil,islope)
3954     &         * mmol(igcm_h2o_vap)/(mugaz*r)
3955             ENDDO
3956          ENDDO
3957        ENDDO
3958
3959      CALL write_output("waterdensity_soil",
3960     &     "rhowater_soil",'kg.m-3',
3961     &     rhowater_soil(:,:,iflat))
3962      CALL write_output("waterdensity_surface",
3963     &     "rhowater_surface",'kg.m-3',
3964     &     rhowater_surf(:,iflat))
3965      DO islope = 1,nslope
3966        write(str2(1:2),'(i2.2)') islope
3967        CALL write_output("waterdensity_soil_slope"//str2,
3968     &     "rhowater_soil_slope"//str2,'kg.m-3',
3969     &     rhowater_soil(:,:,islope))
3970        CALL write_output("waterdensity_surface_slope"//str2,
3971     &     "rhowater_surface"//str2,'kg.m-3',
3972     &     rhowater_surf(:,islope))
3973      ENDDO
3974
3975      CALL write_output("h2o_layer1","h2o mass mr in the first layer",
3976     &     'kg/kg',zq(:,1,igcm_h2o_vap))
3977      CALL write_output("co2_layer1","co2 mass mr in the first layer",
3978     &     'kg/kg',zq(:,1,igcm_co2))
3979        ENDIF ! of IF (water)
3980
3981!PREVIOUSLY IN 1D ONLY
3982
3983c        ==========================================================
3984c        END OF WRITEDIAGFI
3985c        ==========================================================
3986#endif
3987! of ifdef MESOSCALE
3988
3989c      ELSE     ! if(ngrid.eq.1)
3990
3991c#ifndef MESOSCALE
3992c         write(*,
3993c     &    '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
3994c     &    zls*180./pi,odpref,tau_pref_scenario
3995c#endif
3996
3997c      END IF       ! if(ngrid.ne.1)
3998
3999! test for co2 conservation with co2 microphysics
4000      if (igcm_co2_ice.ne.0) then
4001        co2totB = 0. ! added by C.M.
4002        do ig=1,ngrid
4003          do l=1,nlayer
4004            co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
4005     &             (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
4006     &        +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
4007          enddo
4008          co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
4009        enddo
4010      else
4011        co2totB = 0. ! added by C.M.
4012        do ig=1,ngrid
4013          do l=1,nlayer
4014            co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
4015     &             (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
4016          enddo
4017          co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
4018        enddo
4019      endif ! of if (igcm_co2_ice.ne.0)
4020      co2conservation = (co2totA-co2totB)/co2totA
4021        call write_output( 'co2conservation',
4022     &                     'Total CO2 mass conservation in physic',
4023     &                     'kg', co2conservation)
4024! XIOS outputs
4025#ifdef CPP_XIOS
4026      ! Send fields to XIOS: (NB these fields must also be defined as
4027      ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
4028
4029      CALL send_xios_field("controle",tab_cntrl_mod,1)
4030
4031      CALL send_xios_field("ap",ap,1)
4032      CALL send_xios_field("bp",bp,1)
4033      CALL send_xios_field("aps",aps,1)
4034      CALL send_xios_field("bps",bps,1)
4035
4036      if (lastcall.and.is_omp_master) then
4037        write(*,*) "physiq lastcall: call xios_context_finalize"
4038        call xios_context_finalize
4039      endif
4040#endif
4041
4042      if (check_physics_outputs) then
4043        ! Check the validity of updated fields at the end of the physics step
4044        call check_physics_fields("end of physiq:",zt,zu,zv,zplev,zq)
4045      endif
4046
4047      icount=icount+1
4048
4049      END SUBROUTINE physiq
4050
4051      END MODULE physiq_mod
Note: See TracBrowser for help on using the repository browser.