source: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 @ 3661

Last change on this file since 3661 was 3661, checked in by emoisan, 4 months ago

Titan CRM:
Add Titan interface in INTERFACES_V4
Adapt module_model_constants.F to Titan
Add new tracer_mode for Titan (CH4 scalar)
Add new communication of variables between LMDZ.TITAN and WRF
Allow microphysics for Mesoscale in physiq_mod.F90
EMo

File size: 86.5 KB
Line 
1      module physiq_mod
2     
3      implicit none
4     
5      contains
6     
7      subroutine physiq(ngrid,nlayer,nq,   &
8                  nametrac,                &
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 radinc_h, only : L_NSPECTI,L_NSPECTV
17      use radcommon_h, only: sigma, gzlat, grav, BWNV, WAVEI, WAVEV
18      use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot, botCH4
19      use comdiurn_h, only: coslat, sinlat, coslon, sinlon
20      use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
21      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
22      use datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file, nudging_file
23      use geometry_mod, only: latitude, latitude_deg, longitude, cell_area
24      USE comgeomfi_h, only: totarea, totarea_planet
25      USE tracer_h
26      use time_phylmdz_mod, only: ecritphy, iphysiq, nday
27      use phyetat0_mod, only: phyetat0
28      use phyredem, only: physdem0, physdem1
29      use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
30      use mod_phys_lmdz_para, only : is_master
31      use planete_mod, only: apoastr, periastr, year_day, peri_day, &
32                            obliquit, nres, z0
33      use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
34      use time_phylmdz_mod, only: daysec
35#ifndef MESOSCALE
36      use logic_mod, only: moyzon_ch
37      use moyzon_mod, only:  zphibar, zphisbar, zplevbar, zplaybar, &
38                             zzlevbar, zzlaybar, ztfibar, zqfibar
39#endif
40      use callkeys_mod
41      use phys_state_var_mod
42      use turb_mod, only : q2,sensibFlux,turb_resolved
43#ifndef MESOSCALE
44      use vertical_layers_mod, only: presnivs, pseudoalt
45      use mod_phys_lmdz_omp_data, ONLY: is_omp_master
46#else
47use comm_wrf, only : comm_HR_SW, comm_HR_LW, &
48                     comm_FLUXTOP_DN, comm_FLUXABS_SW,&
49                     comm_FLUXTOP_LW, comm_FLUXSURF_SW,&
50                     comm_FLUXSURF_LW, comm_FLXGRD,&
51                     comm_zqfi_omp, comm_zdtlc
52#endif
53#ifdef CPP_XIOS     
54      use xios_output_mod, only: initialize_xios_output, &
55                                 update_xios_timestep, &
56                                 send_xios_field
57      use wxios, only: wxios_context_init, xios_context_finalize
58#endif
59      use muphy_diag
60      implicit none
61
62
63!==================================================================
64!     
65!     Purpose
66!     -------
67!     Central subroutine for all the physics parameterisations in the
68!     universal model. Originally adapted from the Mars LMDZ model.
69!
70!     The model can be run without or with tracer transport
71!     depending on the value of "tracer" in file "callphys.def".
72!
73!
74!   It includes:
75!
76!      I. Initialization :
77!         I.1 Firstcall initializations.
78!         I.2 Initialization for every call to physiq.
79!
80!      II. Compute radiative transfer tendencies (longwave and shortwave) :
81!         II.a Option 1 : Call correlated-k radiative transfer scheme.
82!         II.b Option 2 : Call Newtonian cooling scheme.
83!         II.c Option 3 : Atmosphere has no radiative effect.
84!
85!      III. Vertical diffusion (turbulent mixing) :
86!
87!      IV. Dry Convective adjustment :
88!
89!      V. Tracers
90!         V.1. Microphysics
91!         V.2. Chemistry
92!         V.3. Updates (pressure variations, surface budget).
93!         V.4. Surface Tracer Update.
94!
95!      VI. Surface and sub-surface soil temperature.
96!
97!      VII. Perform diagnostics and write output files.
98!
99!
100!   arguments
101!   ---------
102!
103!   INPUT
104!   -----
105!
106!    ngrid                 Size of the horizontal grid.
107!    nlayer                Number of vertical layers.
108!    nq                    Number of advected fields.
109!    nametrac              Name of corresponding advected fields.
110!
111!    firstcall             True at the first call.
112!    lastcall              True at the last call.
113!
114!    pday                  Number of days counted from the North. Spring equinoxe.
115!    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
116!    ptimestep             timestep (s).
117!
118!    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa).
119!    pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
120!    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2.s-2).
121!
122!    pu(ngrid,nlayer)      u, zonal component of the wind (ms-1).
123!    pv(ngrid,nlayer)      v, meridional component of the wind (ms-1).
124!
125!    pt(ngrid,nlayer)      Temperature (K).
126!
127!    pq(ngrid,nlayer,nq)   Advected fields.
128!
129!    pudyn(ngrid,nlayer)    \
130!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
131!    ptdyn(ngrid,nlayer)     / corresponding variables.
132!    pqdyn(ngrid,nlayer,nq) /
133!    flxw(ngrid,nlayer)      vertical mass flux (kg/s) at layer lower boundary
134!
135!   OUTPUT
136!   ------
137!
138!    pdu(ngrid,nlayer)        \
139!    pdv(ngrid,nlayer)         \  Temporal derivative of the corresponding
140!    pdt(ngrid,nlayer)         /  variables due to physical processes.
141!    pdq(ngrid,nlayer,nq)     /
142!    pdpsrf(ngrid)           /
143!
144!
145!     Authors
146!     -------
147!           Frederic Hourdin        15/10/93
148!           Francois Forget        1994
149!           Christophe Hourdin        02/1997
150!           Subroutine completely rewritten by F. Forget (01/2000)
151!           Water ice clouds: Franck Montmessin (update 06/2003)
152!           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
153!           New correlated-k radiative scheme: R. Wordsworth (2009)
154!           Many specifically Martian subroutines removed: R. Wordsworth (2009)       
155!           Improved water cycle: R. Wordsworth / B. Charnay (2010)
156!           To F90: R. Wordsworth (2010)
157!           New turbulent diffusion scheme: J. Leconte (2012)
158!           Loops converted to F90 matrix format: J. Leconte (2012)
159!           No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
160!           Purge of the code : M. Turbet (2015)
161!           Fork for Titan : J. Vatant d'Ollone (2017)
162!                + clean of all too-generic (ocean, water, co2 ...) routines
163!                + Titan's chemistry
164!           Microphysical moment model: J.Burgalat / B. de Batz de Trenquelléon (2022-2023)
165!           Optics for haze and clouds: B. de Batz de Trenquelléon (2023)
166!============================================================================================
167
168! ---------------
169! 0. DECLARATIONS
170! ---------------
171
172include "netcdf.inc"
173
174! Arguments :
175! -----------
176
177!   INPUTS:
178!   -------
179
180      integer,intent(in) :: ngrid             ! Number of atmospheric columns.
181      integer,intent(in) :: nlayer            ! Number of atmospheric layers.
182      integer,intent(in) :: nq                ! Number of tracers.
183      character*30,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics.
184     
185      logical,intent(in) :: firstcall ! Signals first call to physics.
186      logical,intent(in) :: lastcall  ! Signals last call to physics.
187     
188      real,intent(in) :: pday                  ! Number of elapsed sols since reference Ls=0.
189      real,intent(in) :: ptime                 ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
190      real,intent(in) :: ptimestep             ! Physics timestep (s).
191      real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
192      real,intent(in) :: pplay(ngrid,nlayer)   ! Mid-layer pressure (Pa).
193      real,intent(in) :: pphi(ngrid,nlayer)    ! Geopotential at mid-layer (m2s-2).
194      real,intent(in) :: pu(ngrid,nlayer)      ! Zonal wind component (m/s).
195      real,intent(in) :: pv(ngrid,nlayer)      ! Meridional wind component (m/s).
196      real,intent(in) :: pt(ngrid,nlayer)      ! Temperature (K).
197      real,intent(in) :: pq(ngrid,nlayer,nq)   ! Tracers (kg/kg_of_air).
198      real,intent(in) :: flxw(ngrid,nlayer)    ! Vertical mass flux (ks/s) at lower boundary of layer
199
200!   OUTPUTS:
201!   --------
202
203!     Physical tendencies :
204
205      real,intent(out) :: pdu(ngrid,nlayer)    ! Zonal wind tendencies (m/s/s).
206      real,intent(out) :: pdv(ngrid,nlayer)    ! Meridional wind tendencies (m/s/s).
207      real,intent(out) :: pdt(ngrid,nlayer)    ! Temperature tendencies (K/s).
208      real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
209      real,intent(out) :: pdpsrf(ngrid)        ! Surface pressure tendency (Pa/s).
210
211! Local saved variables:
212! ----------------------
213      integer,save :: day_ini                                      ! Initial date of the run (sol since Ls=0).
214      integer,save :: icount                                       ! Counter of calls to physiq during the run.
215!$OMP THREADPRIVATE(day_ini,icount)
216
217! Local variables :
218! -----------------
219
220      real zh(ngrid,nlayer)               ! Potential temperature (K).
221      real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
222
223      integer l,ig,ierr,iq,nw,isoil,ilat,lat_idx,i,j
224
225      ! FOR DIAGNOSTIC :
226
227      real zls                       ! Solar longitude (radians).
228      real zlss                      ! Sub solar point longitude (radians).
229      real zday                      ! Date (time since Ls=0, calculated in sols).
230      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers (ref : local surf).
231      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries (ref : local surf).
232      real zzlay_eff(ngrid,nlayer)   ! Effective altitude at the middle of the atmospheric layers (ref : geoid ).
233      real zzlev_eff(ngrid,nlayer+1) ! Effective altitude at the atmospheric layer boundaries ( ref : geoid ).
234     
235! TENDENCIES due to various processes :
236
237      ! For Surface Temperature : (K/s)
238      real zdtsurf(ngrid)     ! Cumulated tendencies.
239      real zdtsurfmr(ngrid)   ! Mass_redistribution routine.
240      real zdtsurfevap(ngrid) ! Evaporation.
241      real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
242           
243      ! For Atmospheric Temperatures : (K/s)   
244      real zdtdif(ngrid,nlayer)                       ! Turbdiff/vdifc routines.
245      real zdtmr(ngrid,nlayer)                        ! Mass_redistribution routine.
246      real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
247      real zdtlc(ngrid,nlayer)                        ! Condensation heating rate.
248                             
249      ! For Surface Tracers : (kg/m2/s)
250      real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
251      real zdqsdif(ngrid,nq)                ! Turbdiff/vdifc routines.
252      real zdqsurfmr(ngrid,nq)              ! Mass_redistribution routine.
253                 
254      ! For Tracers : (kg/kg_of_air/s)
255      real zdqadj(ngrid,nlayer,nq)    ! Convadj routine.
256      real zdqdif(ngrid,nlayer,nq)    ! Turbdiff/vdifc routines.
257      real zdqevap(ngrid,nlayer)      ! Turbdiff routine.
258      real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
259     
260      real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency (chemistry routine).
261     
262      real zdqmufi(ngrid,nlayer,nq)   ! Microphysical tendency.
263     
264      real zdqfibar(ngrid,nlayer,nq)   ! For 2D chemistry
265      real zdqmufibar(ngrid,nlayer,nq) ! For 2D chemistry
266                       
267      ! For Winds : (m/s/s)
268      real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine.
269      real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer)   ! Mass_redistribution routine.
270      real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
271      real zdhdif(ngrid,nlayer)                      ! Turbdiff/vdifc routines.
272      real zdhadj(ngrid,nlayer)                      ! Convadj routine.
273      real zdundg(ngrid,nlayer)                      ! Nudging for zonal wind.
274     
275      ! For Pressure and Mass :
276      real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
277      real zdmassmr_col(ngrid)    ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
278      real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
279     
280   
281! Local variables for LOCAL CALCULATIONS:
282! ---------------------------------------
283      real zflubid(ngrid)
284      real zplanck(ngrid),zpopsk(ngrid,nlayer)
285      real ztim1,ztim2,ztim3, z1,z2
286      real ztime_fin
287      real zdh(ngrid,nlayer)
288      real gmplanet
289      real taux(ngrid),tauy(ngrid)
290      real factlat
291      real zundg(nlayer)
292
293
294! local variables for DIAGNOSTICS : (diagfi & stat)
295! -------------------------------------------------
296      real ps(ngrid)                                     ! Surface Pressure.
297      real zt(ngrid,nlayer)                              ! Atmospheric Temperature.
298      real zu(ngrid,nlayer),zv(ngrid,nlayer)             ! Zonal and Meridional Winds.
299      real zq(ngrid,nlayer,nq)                           ! Atmospheric Tracers.
300      real zdtadj(ngrid,nlayer)                          ! Convadj Diagnostic.
301      real zdtdyn(ngrid,nlayer)                          ! Dynamical Heating (K/s).
302      real zdudyn(ngrid,nlayer)                          ! Dynamical Zonal Wind tendency (m.s-2).
303
304      real zhorizwind(ngrid,nlayer) ! Horizontal Wind (sqrt(u*u+v*v))
305
306      real vmr(ngrid,nlayer)                        ! volume mixing ratio
307      real time_phys
308     
309      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
310     
311      ! To test energy conservation (RW+JL)
312      real mass(ngrid,nlayer),massarea(ngrid,nlayer)
313      real dEtot, dEtots, AtmToSurf_TurbFlux
314      real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
315!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
316      real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
317      real dEdiffs(ngrid),dEdiff(ngrid)
318     
319      ! JL12 : conservation test for mean flow kinetic energy has been disabled temporarily
320      real dItot, dItot_tmp, dVtot, dVtot_tmp
321      real dWtot, dWtot_tmp, dWtots, dWtots_tmp
322     
323      ! For Clear Sky Case.
324      real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid)  ! For SW/LW flux.
325      real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid)                            ! For SW/LW flux.
326      real albedo_equivalent1(ngrid)                                         ! For Equivalent albedo calculation.
327      real tf, ntf   
328
329      ! Miscellaneous :
330      character(len=10) :: tmp1
331      character(len=10) :: tmp2
332      character*2 :: str2
333
334!#ifndef MESOSCALE
335
336! Local variables for Titan chemistry and microphysics
337! ----------------------------------------------------
338 
339      real,save :: ctimestep ! Chemistry timestep (s)
340!$OMP THREADPRIVATE(ctimestep)
341 
342      ! Chemical tracers in molar fraction [mol/mol]
343      real, dimension(ngrid,nlayer,nkim) :: ychim
344      real, dimension(ngrid,nlayer,nkim) :: ychimbar ! For 2D chemistry
345
346      ! Molar fraction tendencies (chemistry, condensation and evaporation) for tracers [mol/mol/s]
347      real, dimension(ngrid,nlayer,nq)              :: dyccond    ! Condensation rate. NB : for all tracers, as we want to use indx on it.
348      real, dimension(ngrid,nlayer,size(ices_indx)) :: dmuficond  ! Condensation rate from microphysics [kg/kg/s].
349      real, dimension(ngrid,nlayer,nq)              :: dyccondbar ! For 2D chemistry
350      real, dimension(ngrid)                        :: dycevapCH4 ! Surface "pseudo-evaporation" rate for CH4.
351
352      ! Saturation profiles [mol/mol]
353      real, dimension(ngrid,nlayer,nkim) :: ysat
354     
355      ! Temporary fraction of CH4 [mol/mol]
356      real, dimension(ngrid,nlayer) :: tpq_CH4
357
358      real :: i2e(ngrid,nlayer) ! int 2 ext factor (X.kg-1 -> X.m-3)
359
360#ifdef USE_QTEST
361      real,save,dimension(:,:,:), allocatable :: tpq ! Tracers for decoupled microphysical tests ( temporary in 01/18 )
362!$OMP THREADPRIVATE(tpq)
363      real,dimension(ngrid,nlayer,nq) :: dtpq ! (temporary in 01/18)
364#endif
365
366#ifdef MESOSCALE
367      LOGICAL, SAVE ::  moyzon_ch ! used for zonal averages in Titan
368      REAL, ALLOCATABLE :: zplevbar(:,:)
369      REAL, ALLOCATABLE :: zplaybar(:,:)
370      REAL, ALLOCATABLE :: ztfibar(:,:)
371      REAL, ALLOCATABLE :: zqfibar(:,:,:)
372      REAL, ALLOCATABLE :: zphibar(:,:)
373      REAL, ALLOCATABLE :: zphisbar(:)
374      REAL, ALLOCATABLE :: zzlevbar(:,:)
375      REAL, ALLOCATABLE :: zzlaybar(:,:)
376!     REAL, DIMENSION(:,:)   :: zplevbar(:,:)
377!     REAL, DIMENSION(:,:)   :: zplaybar(:,:)
378!     REAL, DIMENSION(:,:)   :: ztfibar(:,:)
379!     REAL, DIMENSION(:,:,:) :: zqfibar(:,:,:)
380!     REAL, DIMENSION(:,:)   :: zphibar(:,:)
381!     REAL, DIMENSION(:)     :: zphisbar(:)
382!     REAL, DIMENSION(:,:)   :: zzlevbar(:,:)
383!     REAL, DIMENSION(:,:)   :: zzlaybar(:,:)
384#endif
385
386!-----------------------------------------------------------------------------
387    ! Interface to calmufi
388    !   --> needed in order to pass assumed-shape arrays. Otherwise we must put calmufi in a module
389    !       (to have an explicit interface generated by the compiler).
390    !   Or one can put calmufi in MMP_GCM module (in muphytitan).
391    INTERFACE
392      SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdq)
393        REAL(kind=8), INTENT(IN)                 :: dt    !! Physics timestep (s).
394        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev  !! Pressure levels (Pa).
395        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev  !! Altitude levels (m).
396        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play  !! Pressure layers (Pa).
397        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay  !! Altitude at the center of each layer (m).
398        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d   !! Latitude-Altitude depending gravitational acceleration (m.s-2).
399        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K). 
400        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(X.kg^{-1}}\)).
401        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi !! Tendency from former processes for tracers (\(X.kg^{-1}}\)).
402        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq   !! Microphysical tendency for tracers (\(X.kg^{-1}}\)).
403      END SUBROUTINE calmufi
404    END INTERFACE
405
406!#endif
407     
408      logical file_ok
409
410!==================================================================================================
411
412! -----------------
413! I. INITIALISATION
414! -----------------
415
416! --------------------------------
417! I.1   First Call Initialisation.
418! --------------------------------
419      if (firstcall) then
420
421#ifdef USE_QTEST
422        allocate(tpq(ngrid,nlayer,nq))
423        tpq(:,:,:) = pq(:,:,:)
424#endif
425        ! Initialisation of nmicro as well as tracers names, indexes ...
426        if (ngrid.ne.1) then ! Already done in rcm1d
427           call initracer2(nq,nametrac) ! WARNING JB (27/03/2018): should be wrapped in an OMP SINGLE statement (see module notes)
428        endif
429
430        ! Allocate saved arrays (except for 1D model, where this has already been done)
431#ifndef MESOSCALE
432        if (ngrid>1) call phys_state_var_init(nq)
433#endif
434
435!        Variables set to 0
436!        ~~~~~~~~~~~~~~~~~~
437         dtrad(:,:) = 0.D0
438         zdtlc(:,:) = 0.D0
439         fluxrad(:) = 0.D0
440         zdtsw(:,:) = 0.D0
441         zdtlw(:,:) = 0.D0
442         zpopthi(:,:,:,:) = 0.D0
443         zpopthv(:,:,:,:) = 0.D0
444         zpoptti(:,:,:,:) = 0.D0
445         zpopttv(:,:,:,:) = 0.D0
446
447!        Initialize setup for correlated-k radiative transfer
448!        JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics.
449!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
450
451         if (corrk) then
452         
453           call system('rm -f surf_vals_long.out')
454
455           write( tmp1, '(i3)' ) L_NSPECTI
456           write( tmp2, '(i3)' ) L_NSPECTV
457           banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
458           banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
459
460           call setspi            ! Basic infrared properties.
461           call setspv            ! Basic visible properties.
462           call sugas_corrk       ! Set up gaseous absorption properties.
463         
464           OLR_nu(:,:) = 0.D0
465           OSR_nu(:,:) = 0.D0
466           
467           int_dtaui(:,:,:) = 0.D0
468           int_dtauv(:,:,:) = 0.D0
469         
470!#ifndef MESOSCALE
471           IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
472             haze_opt_file=trim(datadir)//'/optical_tables/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT'
473             inquire(file=trim(haze_opt_file),exist=file_ok)
474             if(.not.file_ok) then
475               write(*,*) 'The file ',TRIM(haze_opt_file),' with the haze optical properties'
476               write(*,*) 'was not found by optci.F90 ! Check in ', TRIM(datadir)
477               write(*,*) 'that you have the one corresponding to the given spectral resolution !!'
478               write(*,*) 'Meanwhile I abort ...'
479               call abort
480              endif
481           ENDIF
482!#endif           
483
484         endif
485
486#ifdef MESOSCALE
487         moyzon_ch = .false. !no zonal mean for mesoscale
488#endif
489
490#ifndef MESOSCALE
491!        Initialize names and timestep for chemistry
492!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
493         if (callchim) then
494
495            if (moyzon_ch .and. ngrid.eq.1) then
496               print *, "moyzon_ch=",moyzon_ch," and ngrid=1"
497               print *, "Please desactivate zonal mean for 1D !"
498               stop
499            endif
500           
501            ! Chemistry timestep
502            ctimestep = ptimestep*REAL(ichim)
503
504         endif
505#endif
506
507!        Initialize microphysics.
508!        ~~~~~~~~~~~~~~~~~~~~~~~~
509
510         IF (callmufi) THEN
511           ! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement.
512           call inimufi(ptimestep)
513
514           ! initialize microphysics diagnostics arrays.
515           call ini_diag_arrays(ngrid,nlayer,nice)
516
517         ENDIF
518
519#ifdef CPP_XIOS
520        ! Initialize XIOS context
521        write(*,*) "physiq: call wxios_context_init"
522        CALL wxios_context_init
523#endif
524
525!        Read 'startfi.nc' file.
526!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
527#ifndef MESOSCALE
528         call phyetat0(startphy_file,                                 &
529                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
530                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,tankCH4)                       
531#else
532         emis(:)=0.0
533         q2(:,:)=0.0
534         qsurf(:,:)=0.0
535         tankCH4(:)=0.0
536         day_ini = pday
537#endif
538
539#ifndef MESOSCALE
540         if (.not.startphy_file) then
541           ! additionnal "academic" initialization of physics
542           if (is_master) write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
543           tsurf(:)=pt(:,1)
544           if (is_master) write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!"
545           do isoil=1,nsoilmx
546             tsoil(:,isoil)=tsurf(:)
547           enddo
548           if (is_master) write(*,*) "Physiq: initializing day_ini to pdat !"
549           day_ini=pday
550         endif
551#endif
552
553         if (pday.ne.day_ini) then
554           write(*,*) "ERROR in physiq.F90:"
555           write(*,*) "bad synchronization between physics and dynamics"
556           write(*,*) "dynamics day: ",pday
557           write(*,*) "physics day:  ",day_ini
558           stop
559         endif
560         write (*,*) 'In physiq day_ini =', day_ini
561
562!        Initialize albedo calculation.
563!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
564         albedo(:,:)=0.D0
565          albedo_bareground(:)=0.D0
566         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground)
567         
568!        Initialize orbital calculation.
569!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
570         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
571
572
573         if(tlocked)then
574            print*,'Planet is tidally locked at resonance n=',nres
575            print*,'Make sure you have the right rotation rate!!!'
576         endif
577
578
579!        Initialize soil.
580!        ~~~~~~~~~~~~~~~~
581         if (callsoil) then
582         
583            call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
584                      ptimestep,tsurf,tsoil,capcal,fluxgrd)
585
586         else ! else of 'callsoil'.
587         
588            print*,'WARNING! Thermal conduction in the soil turned off'
589            capcal(:)=1.e6
590            fluxgrd(:)=intheat
591            print*,'Flux from ground = ',intheat,' W m^-2'
592           
593         endif ! end of 'callsoil'.
594         
595         icount=1
596           
597
598!        Initialize surface history variable.
599!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
600         qsurf_hist(:,:)=qsurf(:,:)
601
602!        Initialize variable for dynamical heating and zonal wind tendency diagnostic
603!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
604         ztprevious(:,:)=pt(:,:)
605         zuprevious(:,:)=pu(:,:)
606
607
608         if(meanOLR)then         
609            call system('rm -f rad_bal.out') ! to record global radiative balance.           
610            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.           
611            call system('rm -f h2o_bal.out') ! to record global hydrological balance.
612         endif
613
614         ! Read NewNudging.dat and initialize the nudging for pu
615         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
616         if (nudging_u) then
617            nudging_file=trim(datadir)//'/NewNudging.dat'
618            inquire(file=trim(nudging_file),exist=file_ok)
619            if(.not.file_ok) then
620               write(*,*) 'ERROR : The file ',TRIM(nudging_file),' was not found by physiq_mod.F90 ! Check in ', TRIM(datadir)
621               write(*,*) 'Meanwhile I abort ...'
622               call abort
623            endif
624           
625            open(88,file=nudging_file,form='formatted')
626            do ilat = 1, 49
627               read(88,*) u_ref(ilat,:)
628            enddo
629            close(88)
630         endif
631
632#ifndef MESOSCALE
633         if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
634            call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
635                         ptimestep,pday+nday,time_phys,cell_area,          &
636                         albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
637         endif
638#endif         
639         
640         ! XIOS outputs
641#ifdef CPP_XIOS
642
643         write(*,*) "physiq: call initialize_xios_output"
644         call initialize_xios_output(pday,ptime,ptimestep,daysec, &
645                                     presnivs,pseudoalt,wavei,wavev)
646#endif
647         write(*,*) "physiq: end of firstcall"
648      endif ! end of 'firstcall'
649
650! ------------------------------------------------------
651! I.2   Initializations done at every physical timestep:
652! ------------------------------------------------------
653
654#ifdef CPP_XIOS     
655      ! update XIOS time/calendar
656      call update_xios_timestep
657#endif     
658
659      ! Initialize various variables
660      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
661
662      pdt(:,:)       = 0.D0     
663      zdtsurf(:)     = 0.D0
664      zdtsurfevap(:) = 0.D0
665      pdq(:,:,:)     = 0.D0
666      dqsurf(:,:)    = 0.D0
667      pdu(:,:)       = 0.D0
668      pdv(:,:)       = 0.D0
669      pdpsrf(:)      = 0.D0
670      zflubid(:)     = 0.D0     
671      taux(:)        = 0.D0
672      tauy(:)        = 0.D0
673
674      zday=pday+ptime ! Compute time, in sols (and fraction thereof).
675
676      ! Compute Stellar Longitude (Ls), and orbital parameters.
677      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
678      if (season) then
679         call stellarlong(zday,zls)
680      else
681         call stellarlong(noseason_day,zls)
682      end if
683
684      call orbite(zls,dist_star,declin,right_ascen)
685     
686      if (tlocked) then
687              zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi)
688      elseif (diurnal) then
689              zlss=-2.*pi*(zday-.5)
690      else if(diurnal .eqv. .false.) then
691              zlss=9999.
692      endif
693
694
695      ! Compute variations of g with latitude (oblate case).
696      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
697      if (oblate .eqv. .false.) then     
698         gzlat(:,:) = g         
699      else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then     
700         print*,'I need values for flatten, J2, Rmean and MassPlanet to compute gzlat (else set oblate=.false.)'
701         call abort
702      else
703         gmplanet = MassPlanet*grav*1e24
704         do ig=1,ngrid
705            gzlat(ig,:)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - latitude(ig))))
706         end do
707      endif
708     
709      ! Compute altitudes with the geopotential coming from the dynamics.
710      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
711
712      if (eff_gz .eqv. .false.) then
713     
714        do l=1,nlayer         
715           zzlay(:,l) = pphi(:,l) / gzlat(:,l) ! Reference = local surface
716        enddo
717       
718      else ! In this case we consider variations of g with altitude
719     
720        do l=1,nlayer
721           zzlay(:,l) = g*rad*rad / ( g*rad - ( pphi(:,l) + phisfi(:) )) - rad
722           gzlat(:,l) = g*rad*rad / ( rad + zzlay(:,l) )**2
723        end do
724     
725      endif ! if eff_gz
726
727      zzlev(:,1)=0.
728      zzlev(:,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
729      ! JVO 19 : This altitude is indeed dummy for the GCM and fits ptop=0
730      ! but for upper chemistry that's a pb -> we anyway redefine it just after ..
731
732      do l=2,nlayer
733         do ig=1,ngrid
734            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
735            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
736            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
737         enddo
738      enddo     
739
740      ! Effective altitudes ( eg needed for chemistry ) with correct g, and with reference to the geoid
741      ! JVO 19 : We shall always have correct altitudes in chemistry no matter what's in physics
742      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
743!#ifndef MESOSCALE
744      if (moyzon_ch) then ! Zonal averages
745         
746         zzlaybar(1,:)=g*rad*rad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad   ! reference = geoid
747         zzlevbar(1,1)=zphisbar(1)/g       
748         DO l=2,nlayer
749            z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplaybar(1,l-1)-zplevbar(1,l))
750            z2=(zplevbar(1,l)  +zplaybar(1,l))/(zplevbar(1,l)  -zplaybar(1,l))
751            zzlevbar(1,l)=(z1*zzlaybar(1,l-1)+z2*zzlaybar(1,l))/(z1+z2)
752         ENDDO
753         zzlevbar(1,nlayer+1)=zzlaybar(1,nlayer)+(zzlaybar(1,nlayer)-zzlevbar(1,nlayer))
754
755         DO ig=2,ngrid
756            if (latitude(ig).ne.latitude(ig-1)) then       
757               DO l=1,nlayer   
758                  zzlaybar(ig,l)=g*rad*rad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
759               ENDDO
760               zzlevbar(ig,1)=zphisbar(ig)/g
761               DO l=2,nlayer
762                  z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplaybar(ig,l-1)-zplevbar(ig,l))
763                  z2=(zplevbar(ig,l)  +zplaybar(ig,l))/(zplevbar(ig,l)  -zplaybar(ig,l))
764                  zzlevbar(ig,l)=(z1*zzlaybar(ig,l-1)+z2*zzlaybar(ig,l))/(z1+z2)
765               ENDDO
766               zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))         
767            else
768               zzlaybar(ig,:)=zzlaybar(ig-1,:)
769               zzlevbar(ig,:)=zzlevbar(ig-1,:)
770            endif         
771         ENDDO
772
773      else !  if not moyzon
774!#endif
775     
776        DO ig=1,ngrid
777          DO l=1,nlayer   
778            zzlay_eff(ig,l)=g*rad*rad/(g*rad-(pphi(ig,l)+phisfi(ig)))-rad ! reference = geoid
779          ENDDO
780          zzlev_eff(ig,1)=phisfi(ig)/g
781          DO l=2,nlayer
782            z1=(pplay(ig,l-1)+pplev(ig,l))/ (pplay(ig,l-1)-pplev(ig,l))
783            z2=(pplev(ig,l)  +pplay(ig,l))/(pplev(ig,l)  -pplay(ig,l))
784            zzlev_eff(ig,l)=(z1*zzlay_eff(ig,l-1)+z2*zzlay_eff(ig,l))/(z1+z2)
785          ENDDO
786          zzlev_eff(ig,nlayer+1)=zzlay_eff(ig,nlayer)+(zzlay_eff(ig,nlayer)-zzlev_eff(ig,nlayer))
787        ENDDO
788
789!#ifndef MESOSCALE
790      endif  ! moyzon
791!#endif
792
793      ! -------------------------------------------------------------------------------------
794      ! Compute potential temperature
795      ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
796      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
797      do l=1,nlayer         
798         do ig=1,ngrid
799            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
800            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
801            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/gzlat(ig,l)
802            massarea(ig,l)=mass(ig,l)*cell_area(ig)
803         enddo
804      enddo
805
806     ! Compute vertical velocity (m/s) from vertical mass flux
807     ! w = F / (rho*area) and rho = P/(r*T)
808     ! But first linearly interpolate mass flux to mid-layers
809      do l=1,nlayer-1
810         pw(:,l)=0.5*(flxw(:,l)+flxw(:,l+1))
811      enddo
812      pw(:,nlayer)=0.5*flxw(:,nlayer) ! since flxw(nlayer+1)=0
813      do l=1,nlayer
814         pw(:,l)=(pw(:,l)*r*pt(:,l)) / (pplay(:,l)*cell_area(:))
815      enddo
816
817!---------------------------------
818! II. Compute radiative tendencies
819!---------------------------------
820
821      if (callrad) then
822         if( mod(icount-1,iradia).eq.0.or.lastcall) then
823
824          ! Compute local stellar zenith angles
825          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
826            if (tlocked) then
827            ! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
828               ztim1=SIN(declin)
829               ztim2=COS(declin)*COS(zlss)
830               ztim3=COS(declin)*SIN(zlss)
831
832               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
833                            ztim1,ztim2,ztim3,mu0,fract, flatten)
834
835            elseif (diurnal) then
836               ztim1=SIN(declin)
837               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
838               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
839
840               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
841                            ztim1,ztim2,ztim3,mu0,fract, flatten)
842           
843            else if(diurnal .eqv. .false.) then
844
845               call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten)
846               ! WARNING: this function appears not to work in 1D
847
848            endif
849           
850            ! Eclipse incoming sunlight (e.g. Saturn ring shadowing).       
851            if(rings_shadow) then
852                call call_rings(ngrid, ptime, pday, diurnal)
853            endif   
854
855
856            if (corrk) then
857           
858! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
859! II.a Call correlated-k radiative transfer scheme
860! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
861
862               call call_profilgases(nlayer)
863
864               ! standard callcorrk
865               call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,                      &
866                              albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev,&
867                              pt,tsurf,fract,dist_star,                           &
868                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
869                              fluxsurfabs_sw,fluxtop_lw,                          &
870                              fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,                &
871                              int_dtaui,int_dtauv,zpopthi,zpopthv,zpoptti,zpopttv,&
872                              lastcall)
873
874               ! Radiative flux from the sky absorbed by the surface (W.m-2).
875               GSR=0.0
876               fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)+fluxsurfabs_sw(:)
877
878                            !if(noradsurf)then ! no lower surface; SW flux just disappears
879                            !   GSR = SUM(fluxsurf_sw(:)*cell_area(:))/totarea
880                            !   fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)
881                            !   print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
882                            !endif
883
884               ! Net atmospheric radiative heating rate (K.s-1)
885               dtrad(:,:)=zdtsw(:,:)+zdtlw(:,:)
886               
887            elseif(newtonian)then
888           
889! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
890! II.b Call Newtonian cooling scheme
891! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
892               call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
893
894               zdtsurf(:) = +(pt(:,1)-tsurf(:))/ptimestep
895               ! e.g. surface becomes proxy for 1st atmospheric layer ?
896
897            else
898           
899! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
900! II.c Atmosphere has no radiative effect
901! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
902               fluxtop_dn(:)  = fract(:)*mu0(:)*Fat1AU/dist_star**2
903               if(ngrid.eq.1)then ! / by 4 globally in 1D case...
904                  fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
905               endif
906               fluxsurf_sw(:) = fluxtop_dn(:)
907               print*,'------------WARNING---WARNING------------' ! by MT2015.
908               print*,'You are in corrk=false mode, '
909               print*,'and the surface albedo is taken equal to the first visible spectral value'               
910               
911               fluxsurfabs_sw(:) = fluxtop_dn(:)*(1.-albedo(:,1))
912               fluxrad_sky(:)    = fluxsurfabs_sw(:)
913               fluxtop_lw(:)     = emis(:)*sigma*tsurf(:)**4
914
915               dtrad(:,:)=0.D0 ! no atmospheric radiative heating
916
917            endif ! end of corrk
918
919         endif ! of if(mod(icount-1,iradia).eq.0)
920       
921
922         ! Transformation of the radiative tendencies
923         ! ------------------------------------------
924         zplanck(:)=tsurf(:)*tsurf(:)
925         zplanck(:)=emis(:)*sigma*zplanck(:)*zplanck(:)
926         fluxrad(:)=fluxrad_sky(:)-zplanck(:)
927         pdt(:,:)=pdt(:,:)+dtrad(:,:)
928         
929         ! Test of energy conservation
930         !----------------------------
931         if(enertest)then
932            call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
933            call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW)
934            !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
935            call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
936            call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*cell_area(:)/totarea_planet,dEtotsLW)
937            dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
938            dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
939            if (is_master) then
940               print*,'---------------------------------------------------------------'
941               print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
942               print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
943               print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
944               print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
945               print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
946               print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
947            endif
948         endif ! end of 'enertest'
949
950      endif ! of if (callrad)
951
952
953
954!  --------------------------------------------
955!  III. Vertical diffusion (turbulent mixing) :
956!  --------------------------------------------
957      if (calldifv) then
958     
959         zflubid(:)=fluxrad(:)+fluxgrd(:)
960
961         ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
962         if (UseTurbDiff) then
963         
964            call turbdiff(ngrid,nlayer,nq,                       &
965                          ptimestep,capcal,lwrite,               &
966                          pplay,pplev,zzlay,zzlev,z0,            &
967                          pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
968                          pdt,pdq,zflubid,                       &
969                          zdudif,zdvdif,zdtdif,zdtsdif,          &
970                          sensibFlux,q2,zdqdif,zdqsdif,          &
971                          taux,tauy,lastcall)
972
973         else
974         
975            zdh(:,:)=pdt(:,:)/zpopsk(:,:)
976 
977            call vdifc(ngrid,nlayer,nq,zpopsk,                &
978                       ptimestep,capcal,lwrite,               &
979                       pplay,pplev,zzlay,zzlev,z0,            &
980                       pu,pv,zh,pq,tsurf,emis,qsurf,          &
981                       zdh,pdq,zflubid,                       &
982                       zdudif,zdvdif,zdhdif,zdtsdif,          &
983                       sensibFlux,q2,zdqdif,zdqsdif,          &
984                       taux,tauy,lastcall)
985
986            zdtdif(:,:)=zdhdif(:,:)*zpopsk(:,:) ! for diagnostic only
987            zdqevap(:,:)=0.
988
989         end if !end of 'UseTurbDiff'
990
991         if (.not. turb_resolved) then
992           pdv(:,:)=pdv(:,:)+zdvdif(:,:)
993           pdu(:,:)=pdu(:,:)+zdudif(:,:)
994           pdt(:,:)=pdt(:,:)+zdtdif(:,:)
995         endif
996         
997         zdtsurf(:)=zdtsurf(:)+zdtsdif(:)
998
999         if (tracer) then
1000           pdq(:,:,:)=pdq(:,:,:)+ zdqdif(:,:,:)
1001           dqsurf(:,:)=dqsurf(:,:) + zdqsdif(:,:)
1002         end if ! of if (tracer)
1003
1004
1005         ! test energy conservation
1006         !-------------------------
1007         if(enertest)then
1008         
1009            dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
1010            do ig = 1, ngrid
1011               dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
1012               dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
1013            enddo
1014           
1015            call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
1016            dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
1017            call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
1018            call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
1019           
1020            if (is_master) then
1021           
1022               if (UseTurbDiff) then
1023                         print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
1024                         print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
1025                  print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1026               else
1027                         print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
1028                         print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
1029                         print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
1030               end if
1031            endif ! end of 'is_master'
1032           
1033         ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
1034         endif ! end of 'enertest'
1035
1036      else ! calldifv
1037
1038         if(.not.newtonian)then
1039
1040            zdtsurf(:) = zdtsurf(:) + (fluxrad(:) + fluxgrd(:))/capcal(:)
1041
1042         endif
1043
1044      endif ! end of 'calldifv'
1045
1046
1047!----------------------------------
1048!   IV. Dry convective adjustment :
1049!----------------------------------
1050
1051      if(calladj) then
1052
1053         zdh(:,:) = pdt(:,:)/zpopsk(:,:)
1054         zduadj(:,:)=0.D0
1055         zdvadj(:,:)=0.D0
1056         zdhadj(:,:)=0.D0
1057
1058
1059         call convadj(ngrid,nlayer,nq,ptimestep,            &
1060                      pplay,pplev,zpopsk,                   &
1061                      pu,pv,zh,pq,                          &
1062                      pdu,pdv,zdh,pdq,                      &
1063                      zduadj,zdvadj,zdhadj,                 &
1064                      zdqadj)
1065
1066         pdu(:,:) = pdu(:,:) + zduadj(:,:)
1067         pdv(:,:) = pdv(:,:) + zdvadj(:,:)
1068         pdt(:,:)    = pdt(:,:) + zdhadj(:,:)*zpopsk(:,:)
1069         zdtadj(:,:) = zdhadj(:,:)*zpopsk(:,:) ! for diagnostic only
1070
1071         if(tracer) then
1072            pdq(:,:,:) = pdq(:,:,:) + zdqadj(:,:,:)
1073         end if
1074
1075         ! Test energy conservation
1076         if(enertest)then
1077            call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot)
1078            if (is_master) print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
1079         endif
1080
1081         
1082      endif ! end of 'calladj'
1083     
1084
1085!---------------------------------------------
1086!   V. Specific parameterizations for tracers
1087!---------------------------------------------
1088
1089      if (tracer) then
1090
1091   ! -------------------
1092   !   V.1 Microphysics
1093   ! -------------------
1094
1095         ! We must call microphysics before chemistry, for condensation !
1096         if (callmufi) then
1097            zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on (could be done cleaner)
1098
1099#ifdef USE_QTEST
1100            dtpq(:,:,:) = 0.D0 ! we want tpq to go only through mufi
1101            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
1102            tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
1103
1104#else   
1105            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
1106            pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
1107
1108            ! Sanity check (way safer to be done here rather than within YAMMS)
1109            ! Important : the sanity check intentionally include the former processes tendency !
1110            ! NB : Despite this sanity check there might be still some unphysical values going through :
1111            !        - Negatives, but harmless as it will be only for the output files
1112            !          just remove them in post-proc.
1113            !        - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if
1114            !          you want to compute optics from radii.
1115            WHERE ((pq(:,:,1)+pdq(:,:,1)*ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)*ptimestep < 0.D0))
1116               pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep
1117               pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep
1118            ENDWHERE
1119            WHERE ((pq(:,:,3)+pdq(:,:,3)*ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)*ptimestep < 0.D0))
1120               pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep
1121               pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
1122            ENDWHERE
1123            IF (callclouds) THEN
1124               WHERE ((pq(:,:,5)+pdq(:,:,5)*ptimestep < 0.D0) .OR. (pq(:,:,6)+pdq(:,:,6)*ptimestep < 0.D0))
1125                  pdq(:,:,5) = (epsilon(1.0)-1.D0)*pq(:,:,5)/ptimestep
1126                  pdq(:,:,6) = (epsilon(1.0)-1.D0)*pq(:,:,6)/ptimestep
1127               ENDWHERE
1128               DO iq = 1, size(ices_indx)
1129                  ! For ices :
1130                  WHERE (pq(:,:,ices_indx(iq))+pdq(:,:,ices_indx(iq))*ptimestep < 0.D0)
1131                     pdq(:,:,ices_indx(iq)) = (epsilon(1.0)-1.D0)*pq(:,:,ices_indx(iq))/ptimestep
1132                  ENDWHERE
1133                  ! For gases :
1134                  WHERE (pq(:,:,gazs_indx(iq))+pdq(:,:,gazs_indx(iq))*ptimestep < 0.D0)
1135                     pdq(:,:,gazs_indx(iq)) = (epsilon(1.0)-1.D0)*pq(:,:,gazs_indx(iq))/ptimestep
1136                  ENDWHERE
1137               ENDDO
1138            ENDIF
1139            ! In case there is no clouds, in the troposphere the moments are fixed to evacuate all aerosols
1140            IF (.NOT. callclouds) THEN
1141               WHERE (pplay(:,:) .gt. 1000.)
1142                  pdq(:,:,1) = 0.
1143                  pdq(:,:,2) = 0.
1144               ENDWHERE
1145               WHERE (pplay(:,:) .gt. 5000.)
1146                  pdq(:,:,3) = 0.
1147                  pdq(:,:,4) = 0.
1148               ENDWHERE
1149            ENDIF
1150#endif
1151
1152            ! Microphysics condensation for 2D fields to sent non-saturated fields to photochem
1153            if (callclouds .and. moyzon_ch .and. mod(icount-1,ichim).eq.0) then
1154               zdqfibar(:,:,:) = 0.D0 ! We work in zonal average -> forget processes other than condensation
1155               call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
1156                            gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
1157               ! TODO : Add a sanity check here !
1158            endif
1159         
1160            ! Condensation heating rate :
1161            if (callclouds) then
1162               ! Default value -> no condensation [kg/kg_air/s] :
1163               dmuficond(:,:,:) = 0.D0
1164               do iq = 1, size(ices_indx)
1165                  dmuficond(:,:,iq) = zdqmufi(:,:,gazs_indx(iq))
1166               enddo
1167               call cond_muphy(ngrid,nlayer,pt,dmuficond,zdtlc)
1168               !pdt(:,:) = pdt(:,:) + zdtlc(:,:)
1169            endif
1170         endif ! callmufi
1171     
1172#ifndef MESOSCALE
1173  ! -----------------
1174  !   V.2. Chemistry
1175  ! -----------------
1176  !   NB : Must be call last ( brings fields back to an equilibrium )
1177         if (callchim) then
1178
1179            ! o. Convert updated tracers to molar fraction
1180            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1181            do iq = 1,nkim
1182               ychim(:,:,iq) = (pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro)*ptimestep) / rat_mmol(iq+nmicro)
1183            enddo
1184
1185            ! JVO 05/18 : We update zonally averaged fields with condensation
1186            ! as it is compulsory to have correct photochem production. But for other
1187            ! processes ( convadj ... ) we miss them in any case as we work in zonally/diurnal
1188            ! mean -> no fine diurnal/short time evolution, only seasonal evolution only.
1189            if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
1190              do iq = 1,nkim
1191                ychimbar(:,:,iq) =  zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
1192                  if ( callclouds ) then
1193                    ychimbar(:,:,iq) =  ychimbar(:,:,iq) + ( zdqmufibar(:,:,iq+nmicro)*ptimestep / rat_mmol(iq+nmicro) )
1194                  endif
1195              enddo
1196            endif
1197
1198            ! i. Condensation of the 3D tracers after the transport
1199            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1200           
1201            call calc_ysat(ngrid,nlayer,pplay/100.0,pt,ysat) ! Compute saturation profiles for every grid point (!!p in mbar!!)
1202
1203            dyccond(:,:,:) = 0.D0 ! Default value -> no condensation
1204
1205            do iq=1,nkim
1206               where (ychim(:,:,iq).gt.ysat(:,:,iq))
1207                  dyccond(:,:,iq+nmicro) = (-ychim(:,:,iq)+ysat(:,:,iq)) / ptimestep
1208               endwhere
1209            enddo
1210
1211            if (callclouds) then
1212               do iq = 1, size(ices_indx)
1213                  dyccond(:,:,gazs_indx(iq)) = 0.D0 ! Condensation have been calculated in the cloud microphysics
1214               enddo
1215            endif
1216
1217            do iq=1,nkim
1218              ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro)*ptimestep ! update molar ychim for following calchim
1219              pdq(:,:,iq+nmicro) = pdq(:,:,iq+nmicro) + dyccond(:,:,iq+nmicro)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
1220            enddo
1221           
1222
1223            ! ii. 2D zonally averaged fields need to condense and evap before photochem
1224            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1225            if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
1226
1227               call calc_ysat(ngrid,nlayer,zplaybar/100.0,ztfibar,ysat) ! Compute saturation profiles for every grid point for the zon-ave fields
1228
1229               dyccondbar(:,:,:) = 0.D0 ! Default value -> no condensation
1230               
1231               do iq = 1,nkim
1232                     where ( ychimbar(:,:,iq).gt.ysat(:,:,iq) )
1233                        dyccondbar(:,:,iq+nmicro) = ( -ychimbar(:,:,iq)+ysat(:,:,iq) ) / ptimestep
1234                     endwhere
1235               enddo
1236
1237               if (callclouds) then
1238                  do iq = 1, size(ices_indx)
1239                     dyccondbar(:,:,gazs_indx(iq)) = 0.D0 ! Condensation have been calculated in the cloud microphysics
1240                  enddo
1241               endif
1242
1243               do iq=1,nkim
1244                  ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)*ptimestep
1245               enddo
1246
1247            endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq.0 )
1248
1249            ! iii. Photochemistry ( must be call after condensation (and evap of 2D) )
1250            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1251            if( mod(icount-1,ichim).eq.0. ) then
1252               print *, "We enter in the photochemistry ..."
1253
1254               if (moyzon_ch) then ! 2D zonally averaged chemistry
1255                 ! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
1256                 call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,zphisbar,  &
1257                              zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
1258               
1259               else ! 3D chemistry (or 1D run)
1260                 call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,phisfi,  &
1261                              pplay,pplev,zzlay_eff,zzlev_eff,dycchi)
1262               endif ! if moyzon
1263            endif ! if (mod(icount-1,ichim).eq.0)
1264           
1265            ! iv. Surface pseudo-evaporation
1266            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1267            ! Infinite tank of CH4
1268            if (resCH4_inf) then
1269               do ig=1,ngrid
1270                  if ((ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4) then ! + dycchi because ychim not yet updated
1271                     dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
1272                  else
1273                     dycevapCH4(ig) = 0.D0
1274                  endif
1275               enddo
1276
1277            else
1278               ! Fill lakes with precipitation :
1279               !if (REAL(latitude_deg(ig)) .ge. 70.0) then
1280               !   tankCH4(ig) = 200.0 ! [m]
1281               !else if (REAL(latitude_deg(ig)) .le. -70.0) then
1282               !   tankCH4(ig) = 50.0  ! [m]
1283               !else
1284               !   tankCH4(ig) = 0.0   ! [m]
1285               !endif
1286               tankCH4(:) = tankCH4(:) + (mmd_ice_prec(:,1) / 422. * ptimestep) ! [m]
1287               
1288               ! Evaporation of lakes :
1289               if (moyzon_ch) then
1290                 tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
1291               else
1292                 tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep    ! + dycchi because ychim not yet updated [mol/mol]
1293               endif
1294               call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
1295                            pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)
1296               zdtsurf(:) = zdtsurf(:) + zdtsurfevap(:)
1297            endif
1298
1299
1300            pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro) ! convert tendencies to mass mixing ratio
1301           
1302            ! v. Updates and positivity check
1303            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1304            zdqchi(:,:,:)  = 0.D0 ! -> dycchi is saved but for the nmicro tracers we must update to 0 at each step
1305
1306            do iq=1,nkim
1307              zdqchi(:,:,iq+nmicro) = dycchi(:,:,iq)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
1308
1309              where( (pq(:,:,iq+nmicro) + ( pdq(:,:,iq+nmicro)+zdqchi(:,:,iq+nmicro) )*ptimestep ) .LT. 0.)  & ! When using zonal means we set the same tendency
1310                      zdqchi(:,:,iq+nmicro) = 1.D-30 - pdq(:,:,iq+nmicro) - pq(:,:,iq+nmicro)/ptimestep        ! everywhere in longitude -> could lead to negs !
1311            enddo
1312
1313            pdq(:,:,:) = pdq(:,:,:) + zdqchi(:,:,:)
1314
1315         endif ! end of 'callchim'
1316
1317!! END ifndef MESOSCALE
1318#endif
1319
1320  ! ---------------
1321  !   V.3 Updates
1322  ! ---------------
1323
1324         ! Updating Atmospheric Mass and Tracers budgets.
1325         if(mass_redistrib) then
1326
1327            zdmassmr(:,:) = mass(:,:) * zdqevap(:,:)
1328
1329            do ig = 1, ngrid
1330               zdmassmr_col(ig)=SUM(zdmassmr(ig,:))
1331            enddo
1332           
1333            call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
1334            call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
1335            call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
1336
1337            call mass_redistribution(ngrid,nlayer,nq,ptimestep,                     &
1338                                     capcal,pplay,pplev,pt,tsurf,pq,qsurf,          &
1339                                     pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
1340                                     zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
1341         
1342            pdq(:,:,:)  = pdq(:,:,:)  + zdqmr(:,:,:)
1343            dqsurf(:,:) = dqsurf(:,:) + zdqsurfmr(:,:)
1344            pdt(:,:)    = pdt(:,:)    + zdtmr(:,:)
1345            pdu(:,:)    = pdu(:,:)    + zdumr(:,:)
1346            pdv(:,:)    = pdv(:,:)    + zdvmr(:,:)
1347            pdpsrf(:)   = pdpsrf(:)   + zdpsrfmr(:)
1348            zdtsurf(:)  = zdtsurf(:)  + zdtsurfmr(:)
1349         endif
1350
1351  ! -----------------------------
1352  !   V.4. Surface Tracer Update
1353  ! -----------------------------
1354
1355        qsurf(:,:) = qsurf(:,:) + ptimestep*dqsurf(:,:)
1356
1357         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
1358         ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
1359         qsurf_hist(:,:) = qsurf(:,:)
1360
1361      endif! end of if 'tracer'
1362
1363
1364!------------------------------------------------
1365!   VI. Surface and sub-surface soil temperature
1366!------------------------------------------------
1367
1368
1369      ! Increment surface temperature
1370
1371      tsurf(:)=tsurf(:)+ptimestep*zdtsurf(:)
1372
1373      ! Compute soil temperatures and subsurface heat flux.
1374      if (callsoil) then
1375         call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
1376                   ptimestep,tsurf,tsoil,capcal,fluxgrd)     
1377      endif
1378
1379
1380      ! Test energy conservation
1381      if(enertest)then
1382         call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)         
1383         if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
1384      endif
1385
1386
1387!---------------------------------------------------
1388!   VII. Perform diagnostics and write output files
1389!---------------------------------------------------
1390     
1391      ! Nudging of zonal wind !
1392      ! ~~~~~~~~~~~~~~~~~~~~~~~
1393      if (nudging_u) then
1394         zdundg(:,:) = 0.D0
1395         j=1
1396         ! boucle sur les points de grille :
1397         do i = 1, ngrid
1398             ! interpolation linéaire des données dans le fichier lu (sur 49 latitudes)
1399             do while ((u_ref(j,1).ge.latitude_deg(i)).and.(j.lt.49))
1400               j=j+1
1401             enddo
1402             factlat =  (latitude_deg(i)-u_ref(j,1))/(u_ref(j-1,1)-u_ref(j,1))
1403             ! Nudging of the first 23 layers only !!
1404             ! IF CHANGE IN VERTICAL RESOLUTION IN THE FIRST 23 LEVELS, IT DOES NOT WORK !!!
1405             zundg(1:23) = factlat*u_ref(j-1,2:24)+(1-factlat)*u_ref(j,2:24)
1406             zdundg(i,1:23) = (zundg(1:23) - (pu(i,1:23)+pdu(i,1:23)*ptimestep)) / nudging_dt
1407         enddo
1408         
1409         pdu(:,:) = pdu(:,:) + zdundg(:,:)
1410      endif
1411       
1412      ! Note : For output only: the actual model integration is performed in the dynamics.
1413 
1414      ! Temperature, zonal and meridional winds.
1415      zt(:,:) = pt(:,:) + pdt(:,:)*ptimestep
1416      zu(:,:) = pu(:,:) + pdu(:,:)*ptimestep
1417      zv(:,:) = pv(:,:) + pdv(:,:)*ptimestep
1418
1419      ! Diagnostic.
1420      zdtdyn(:,:)     = (pt(:,:)-ztprevious(:,:)) / ptimestep
1421      ztprevious(:,:) = zt(:,:)
1422
1423      zdudyn(:,:)     = (pu(:,:)-zuprevious(:,:)) / ptimestep
1424      zuprevious(:,:) = zu(:,:)
1425
1426      if(firstcall)then
1427         zdtdyn(:,:)=0.D0
1428         zdudyn(:,:)=0.D0
1429      endif
1430
1431      ! Horizotal wind
1432      zhorizwind(:,:) = sqrt( zu(:,:)*zu(:,:) + zv(:,:)*zv(:,:) )
1433
1434      ! Dynamical heating diagnostic.
1435      do ig=1,ngrid
1436         fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp
1437      enddo
1438
1439      ! [Forcage de la photochimie pour les nuages]
1440      if (callclouds) then
1441         do ig = 1, ngrid
1442            do iq = 1, size(ices_indx)
1443               ! C2H2 :
1444               !-------
1445               if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H2") then
1446                  pdq(ig,nlayer-3:,gazs_indx(iq)) = (4.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
1447               endif
1448               ! C2H6 :
1449               !-------
1450               if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H6") then
1451                  pdq(ig,nlayer-3:,gazs_indx(iq)) = (8.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
1452               endif
1453               ! HCN :
1454               !------
1455               if(trim(nameOfTracer(gazs_indx(iq))) .eq. "HCN") then
1456                  pdq(ig,nlayer-3:,gazs_indx(iq)) = (2.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
1457               endif
1458               ! AC6H6 :
1459               !--------
1460               if(trim(nameOfTracer(gazs_indx(iq))) .eq. "AC6H6") then
1461                  pdq(ig,nlayer-3:,gazs_indx(iq)) = (2.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
1462               endif
1463            enddo
1464         enddo
1465      endif
1466
1467      zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
1468
1469      ! Surface pressure.
1470      ps(:) = pplev(:,1) + pdpsrf(:)*ptimestep
1471
1472
1473      ! Surface and soil temperature information
1474      call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1)
1475      call planetwide_minval(tsurf(:),Ts2)
1476      call planetwide_maxval(tsurf(:),Ts3)
1477      if(callsoil)then
1478         TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
1479         if (is_master) then
1480            print*,'          ave[Tsurf]             min[Tsurf]             max[Tsurf]             ave[Tdeep]'
1481            print*,Ts1,Ts2,Ts3,TsS
1482         end if
1483      else
1484         if (is_master) then
1485            print*,'          ave[Tsurf]             min[Tsurf]             max[Tsurf]'
1486            print*,Ts1,Ts2,Ts3
1487         endif
1488      end if
1489
1490
1491      ! Check the energy balance of the simulation during the run
1492      if(corrk)then
1493
1494         call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR)
1495         call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR)
1496         call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR)
1497         call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND)
1498         call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN)
1499         do ig=1,ngrid
1500            if(fluxtop_dn(ig).lt.0.0)then
1501               print*,'fluxtop_dn has gone crazy'
1502               print*,'fluxtop_dn=',fluxtop_dn(ig)
1503               print*,'temp=   ',pt(ig,:)
1504               print*,'pplay=  ',pplay(ig,:)
1505               call abort
1506            endif
1507         end do
1508                     
1509         if(ngrid.eq.1)then
1510            DYN=0.0
1511         endif
1512         
1513         if (is_master) then
1514            print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
1515            print*, ISR,ASR,OLR,GND,DYN
1516         endif
1517
1518         if(enertest .and. is_master)then
1519            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1520            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1521            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1522         endif
1523
1524         if(meanOLR .and. is_master)then
1525            if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
1526               ! to record global radiative balance
1527               open(92,file="rad_bal.out",form='formatted',position='append')
1528               write(92,*) zday,ISR,ASR,OLR
1529               close(92)
1530               open(93,file="tem_bal.out",form='formatted',position='append')
1531               if(callsoil)then
1532                         write(93,*) zday,Ts1,Ts2,Ts3,TsS
1533               else
1534                  write(93,*) zday,Ts1,Ts2,Ts3
1535               endif
1536               close(93)
1537            endif
1538         endif
1539
1540      endif ! end of 'corrk'
1541
1542
1543      ! Diagnostic to test radiative-convective timescales in code.
1544      if(testradtimes)then
1545         open(38,file="tau_phys.out",form='formatted',position='append')
1546         ig=1
1547         do l=1,nlayer
1548            write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1549         enddo
1550         close(38)
1551         print*,'As testradtimes enabled,'
1552         print*,'exiting physics on first call'
1553         call abort
1554      endif
1555
1556      if (is_master) print*,'--> Ls =',zls*180./pi
1557     
1558     
1559!----------------------------------------------------------------------
1560!        Writing NetCDF file  "RESTARTFI" at the end of the run
1561!----------------------------------------------------------------------
1562
1563!        Note: 'restartfi' is stored just before dynamics are stored
1564!              in 'restart'. Between now and the writting of 'restart',
1565!              there will have been the itau=itau+1 instruction and
1566!              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1567!              thus we store for time=time+dtvr
1568
1569
1570
1571      if(lastcall) then
1572         ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1573
1574#ifndef MESOSCALE
1575         if (ngrid.ne.1) then
1576            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1577           
1578            call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
1579                          ptimestep,ztime_fin,                    &
1580                          tsurf,tsoil,emis,q2,qsurf_hist,tankCH4)
1581         endif
1582#endif
1583
1584    endif ! end of 'lastcall'
1585
1586#ifndef MESOSCALE
1587!-----------------------------------------------------------------------------------------------------
1588!           OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
1589!
1590!             Note 1 : output with  period "ecritphy", set in "run.def"
1591!
1592!             Note 2 : writediagfi can also be called from any other subroutine for any variable,
1593!                      but its preferable to keep all the calls in one place ...
1594!-----------------------------------------------------------------------------------------------------
1595
1596
1597      call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
1598      call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
1599      call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
1600      call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
1601      call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1602      call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
1603      call writediagfi(ngrid,"temp","temperature","K",3,zt)
1604      call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
1605      call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
1606      call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
1607      call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
1608      call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
1609     
1610      ! Subsurface temperatures
1611      !call writediagsoil(ngrid,"tempsoil","temperature soil","K",3,tsoil)
1612
1613     ! Total energy balance diagnostics
1614      if(callrad.and.(.not.newtonian))then
1615         call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
1616         call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1617         call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1618         call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1619         
1620         call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
1621         call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
1622         call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
1623         call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
1624         call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
1625         call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
1626
1627         call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)         
1628         call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
1629      endif ! end of 'callrad'
1630       
1631      if(enertest) then
1632         if (calldifv) then
1633            call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
1634            call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
1635            call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
1636            call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
1637            call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
1638         endif
1639         
1640         if (corrk) then
1641            call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
1642            call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
1643         endif
1644      endif ! end of 'enertest'
1645
1646      ! Diagnostics of optical thickness
1647      ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
1648      if (diagdtau) then               
1649         do nw=1,L_NSPECTV
1650            write(str2,'(i2.2)') nw
1651            call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
1652         enddo
1653         do nw=1,L_NSPECTI
1654            write(str2,'(i2.2)') nw
1655            call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
1656         enddo
1657      endif
1658
1659      ! Temporary inclusions for winds diagnostics.
1660      call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
1661      call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
1662
1663      ! Temporary inclusions for heating diagnostics.
1664      call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
1665      call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
1666      call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
1667      call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1668       
1669      ! For Debugging.
1670      call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
1671
1672      ! Output tracers.
1673      if (tracer) then
1674
1675         if (callmufi) then
1676           ! Microphysical tracers are expressed in unit/m3.
1677           ! convert X.kg-1 --> X.m-3 (whereas for optics was -> X.m-2)
1678           i2e(:,:) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer) /(zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
1679
1680#ifdef USE_QTEST
1681            ! Microphysical tracers passed through dyn+phys(except mufi)
1682            call writediagfi(ngrid,"mu_m0as_dp","Dynphys only spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
1683            call writediagfi(ngrid,"mu_m3as_dp","Dynphys only spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
1684            call writediagfi(ngrid,"mu_m0af_dp","Dynphys only fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
1685            call writediagfi(ngrid,"mu_m3af_dp","Dynphys only fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
1686            ! Microphysical tracers passed through mufi only
1687            call writediagfi(ngrid,"mu_m0as_mo","Mufi only spherical mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(1))*i2e)
1688            call writediagfi(ngrid,"mu_m3as_mo","Mufi only spherical mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(2))*i2e)
1689            call writediagfi(ngrid,"mu_m0af_mo","Mufi only fractal mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(3))*i2e)
1690            call writediagfi(ngrid,"mu_m3af_mo","Mufi only fractal mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(4))*i2e)
1691#else
1692            call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
1693            call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
1694            call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
1695            call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
1696#endif
1697           
1698            ! Microphysical diagnostics
1699            call writediagfi(ngrid,"mmd_aer_prec","Total aerosols precipitations",'m',2,mmd_aer_prec)
1700            call writediagfi(ngrid,"mmd_aer_s_flux","Spherical aerosols sedimentation flux",'kg.m-2.s-1',3,mmd_aer_s_flux)
1701            call writediagfi(ngrid,"mmd_aer_f_flux","Fractal aerosols sedimentation flux",'kg.m-2.s-1',3,mmd_aer_f_flux)
1702            call writediagfi(ngrid,"mmd_rc_sph","Spherical mode caracteristic radius",'m',3,mmd_rc_sph)
1703            call writediagfi(ngrid,"mmd_rc_fra","Fractal mode caracteristic radius",'m',3,mmd_rc_fra)
1704         endif ! end of 'callmufi'
1705
1706         ! Chemical tracers
1707         if (callchim) then
1708           do iq=1,nkim
1709             call writediagfi(ngrid,cnames(iq),cnames(iq),'mol/mol',3,zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro))
1710             call writediagfi(ngrid,'dqcond_'//cnames(iq),'dqcond_'//cnames(iq),'mol/mol/s',3,dyccond(:,:,iq+nmicro))
1711           enddo
1712           call writediagfi(ngrid,"evapCH4","Surface CH4 pseudo-evaporation rate",'mol/mol/s',2,dycevapCH4)
1713         endif
1714
1715       endif ! end of 'tracer'
1716
1717#ifdef CPP_XIOS
1718!-----------------------------------------------------------------------------------------------------
1719! XIOS outputs
1720!-----------------------------------------------------------------------------------------------------
1721      ! Send fields to XIOS: (NB these fields must also be defined as
1722      ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
1723
1724      !--------------------------------------------------------
1725      ! General diagnostics :
1726      !--------------------------------------------------------
1727      CALL send_xios_field("ls",zls*180./pi)
1728      CALL send_xios_field("lss",zlss*180./pi)
1729      CALL send_xios_field("RA",right_ascen*180./pi)
1730      CALL send_xios_field("Declin",declin*180./pi)
1731     
1732      ! Atmosphere (3D) :
1733      CALL send_xios_field("temp",zt)
1734      CALL send_xios_field("teta",zh)
1735      CALL send_xios_field("p",pplay)
1736      CALL send_xios_field("u",zu)
1737      CALL send_xios_field("v",zv)
1738      CALL send_xios_field("w",pw)
1739
1740      CALL send_xios_field("area",cell_area)
1741      CALL send_xios_field("pphi",pphi)
1742
1743      ! Surface (2D) :
1744      CALL send_xios_field("ps",ps)
1745      CALL send_xios_field("tsurf",tsurf)
1746      CALL send_xios_field("pphis",phisfi)
1747
1748      ! Total energy balance diagnostics (2D) :
1749      IF (callrad.and.(.not.newtonian)) THEN
1750         CALL send_xios_field("ISR_TOA",fluxtop_dn)
1751         CALL send_xios_field("OLR_TOA",fluxtop_lw)
1752      ENDIF
1753
1754      !--------------------------------------------------------
1755      ! Winds trends :
1756      !--------------------------------------------------------
1757      ! Atmosphere (3D) :
1758      ! du_tot = zdudyn + pdu
1759      CALL send_xios_field("dudyn",zdudyn)
1760      CALL send_xios_field("pdu",pdu)
1761     
1762      ! pdu = zdudif + zduadj + zdundg
1763      CALL send_xios_field("dudif",zdudif)
1764      CALL send_xios_field("duadj",zduadj)
1765      IF (nudging_u) THEN
1766         CALL send_xios_field("dundg",zdundg)
1767      ENDIF
1768     
1769      ! zhorizwind = sqrt(u*u + v*v)
1770      CALL send_xios_field("horizwind",zhorizwind)
1771
1772
1773      !--------------------------------------------------------
1774      ! Heating trends :
1775      !--------------------------------------------------------
1776      ! Atmosphere (3D) :
1777      ! dt_tot = dtdyn + pdt
1778      CALL send_xios_field("dtdyn",zdtdyn)
1779      CALL send_xios_field("pdt",pdt)
1780
1781      ! pdt = dtrad + zdtdif + dtadj + zdtlc
1782      CALL send_xios_field("dtrad",dtrad)
1783      CALL send_xios_field("dtdif",zdtdif)
1784      CALL send_xios_field("dtadj",zdtadj(:,:))
1785      IF (callclouds) THEN
1786         CALL send_xios_field("dtlc",zdtlc)
1787      ENDIF
1788
1789      ! dtrad = zdtsw + zdtlw
1790      CALL send_xios_field("dtsw",zdtsw)
1791      CALL send_xios_field("dtlw",zdtlw)
1792
1793      ! Surface (2D) :
1794      IF(enertest) THEN
1795         IF (calldifv) THEN
1796            CALL send_xios_field("sensibFlux",sensibFlux)
1797         ENDIF
1798      ENDIF
1799      CALL send_xios_field("fluxsurf_lw",fluxsurf_lw(:))
1800      CALL send_xios_field("fluxsurfabs_sw",fluxsurfabs_sw(:))
1801      CALL send_xios_field("emis",emis(:))
1802     
1803      ! dtsurf = dtsdif + dtsurfevap
1804      CALL send_xios_field("dtsurf",zdtsurf(:))
1805      CALL send_xios_field("dtsdif",zdtsdif(:))
1806      CALL send_xios_field("dtsurfevap",zdtsurfevap(:))
1807     
1808
1809      !--------------------------------------------------------
1810      ! Optical diagnostics :
1811      !--------------------------------------------------------
1812      ! Haze opacity :
1813      CALL send_xios_field('ttauhv_08',zpopthv(:,:,8,2))  ! 08 --> 1.983 um
1814      CALL send_xios_field('ttauhv_15',zpopthv(:,:,15,2)) ! 15 --> 1.000 um
1815      CALL send_xios_field('ttauhv_20',zpopthv(:,:,20,2)) ! 20 --> 0.671 um
1816      CALL send_xios_field('ttauhv_23',zpopthv(:,:,23,2)) ! 23 --> 0.346 um
1817      CALL send_xios_field('ttauhi_02',zpopthi(:,:,2,2))  ! 02 --> 175.3 um
1818      CALL send_xios_field('ttauhi_17',zpopthi(:,:,17,2)) ! 17 --> 11.00 um
1819      CALL send_xios_field('ttauhi_23',zpopthi(:,:,23,2)) ! 23 --> 4.849 um
1820      ! Haze extinction :
1821      CALL send_xios_field('kkhv_08',zpopthv(:,:,8,3))
1822      CALL send_xios_field('kkhv_15',zpopthv(:,:,15,3))
1823      CALL send_xios_field('kkhv_20',zpopthv(:,:,20,3))
1824      CALL send_xios_field('kkhv_23',zpopthv(:,:,23,3))
1825      CALL send_xios_field('kkhi_02',zpopthi(:,:,2,3))
1826      CALL send_xios_field('kkhi_17',zpopthi(:,:,17,3))
1827      CALL send_xios_field('kkhi_23',zpopthi(:,:,23,3))
1828      ! Haze single scattering albedo :
1829      CALL send_xios_field('wwhv_08',zpopthv(:,:,8,4))
1830      CALL send_xios_field('wwhv_15',zpopthv(:,:,15,4))
1831      CALL send_xios_field('wwhv_20',zpopthv(:,:,20,4))
1832      CALL send_xios_field('wwhv_23',zpopthv(:,:,23,4))
1833      CALL send_xios_field('wwhi_02',zpopthi(:,:,2,4))
1834      CALL send_xios_field('wwhi_17',zpopthi(:,:,17,4))
1835      CALL send_xios_field('wwhi_23',zpopthi(:,:,23,4))
1836      ! Haze asymmetry parameter :
1837      CALL send_xios_field('gghv_08',zpopthv(:,:,8,5))
1838      CALL send_xios_field('gghv_15',zpopthv(:,:,15,5))
1839      CALL send_xios_field('gghv_20',zpopthv(:,:,20,5))
1840      CALL send_xios_field('gghv_23',zpopthv(:,:,23,5))
1841      CALL send_xios_field('gghi_02',zpopthi(:,:,2,5))
1842      CALL send_xios_field('gghi_17',zpopthi(:,:,17,5))
1843      CALL send_xios_field('gghi_23',zpopthi(:,:,23,5))
1844
1845      ! Diagnostics for haze and clouds :
1846      IF (callclouds) THEN
1847         ! Opacity :
1848         CALL send_xios_field('ttauv_08',zpopttv(:,:,8,2))  ! 08 --> 1.983 um
1849         CALL send_xios_field('ttauv_15',zpopttv(:,:,15,2)) ! 15 --> 1.000 um
1850         CALL send_xios_field('ttauv_20',zpopttv(:,:,20,2)) ! 20 --> 0.671 um
1851         CALL send_xios_field('ttauv_23',zpopttv(:,:,23,2)) ! 23 --> 0.346 um
1852         CALL send_xios_field('ttaui_02',zpoptti(:,:,2,2))  ! 02 --> 175.3 um
1853         CALL send_xios_field('ttaui_17',zpoptti(:,:,17,2)) ! 17 --> 11.00 um
1854         CALL send_xios_field('ttaui_23',zpoptti(:,:,23,2)) ! 23 --> 4.849 um
1855         ! Extinction :
1856         CALL send_xios_field('kkv_08',zpopttv(:,:,8,3))
1857         CALL send_xios_field('kkv_15',zpopttv(:,:,15,3))
1858         CALL send_xios_field('kkv_20',zpopttv(:,:,20,3))
1859         CALL send_xios_field('kkv_23',zpopttv(:,:,23,3))
1860         CALL send_xios_field('kki_02',zpoptti(:,:,2,3))
1861         CALL send_xios_field('kki_17',zpoptti(:,:,17,3))
1862         CALL send_xios_field('kki_23',zpoptti(:,:,23,3))
1863         ! Single scattering albedo :
1864         CALL send_xios_field('wwv_08',zpopttv(:,:,8,4))
1865         CALL send_xios_field('wwv_15',zpopttv(:,:,15,4))
1866         CALL send_xios_field('wwv_20',zpopttv(:,:,20,4))
1867         CALL send_xios_field('wwv_23',zpopttv(:,:,23,4))
1868         CALL send_xios_field('wwi_02',zpoptti(:,:,2,4))
1869         CALL send_xios_field('wwi_17',zpoptti(:,:,17,4))
1870         CALL send_xios_field('wwi_23',zpoptti(:,:,23,4))
1871         ! Asymmetry parameter :
1872         CALL send_xios_field('ggv_08',zpopttv(:,:,8,5))
1873         CALL send_xios_field('ggv_15',zpopttv(:,:,15,5))
1874         CALL send_xios_field('ggv_20',zpopttv(:,:,20,5))
1875         CALL send_xios_field('ggv_23',zpopttv(:,:,23,5))
1876         CALL send_xios_field('ggi_02',zpoptti(:,:,2,5))
1877         CALL send_xios_field('ggi_17',zpoptti(:,:,17,5))
1878         CALL send_xios_field('ggi_23',zpoptti(:,:,23,5))
1879         ! DRAYAER, TAUGAS, DCONT :
1880         CALL send_xios_field('drayaerv_20',zpopttv(:,:,20,6)) ! 20 --> 0.671um
1881         CALL send_xios_field('taugasv_20',zpopttv(:,:,20,7))
1882         CALL send_xios_field('dcontv_20',zpopttv(:,:,20,8))
1883         CALL send_xios_field('drayaeri_17',zpoptti(:,:,17,6)) ! 17 --> 11.00um
1884         CALL send_xios_field('taugasi_17',zpoptti(:,:,17,7))
1885         CALL send_xios_field('dconti_17',zpoptti(:,:,17,8))
1886      ENDIF     
1887
1888      ! Diagnostics for haze and clouds (4D) :
1889      CALL send_xios_field('dtauhi',zpopthi(:,:,:,1))
1890      CALL send_xios_field('tauhi',zpopthi(:,:,:,2))
1891      CALL send_xios_field('khi',zpopthi(:,:,:,3))
1892      CALL send_xios_field('whi',zpopthi(:,:,:,4))
1893      CALL send_xios_field('ghi',zpopthi(:,:,:,5))
1894      CALL send_xios_field('dtauhv',zpopthv(:,:,:,1))
1895      CALL send_xios_field('tauhv',zpopthv(:,:,:,2))
1896      CALL send_xios_field('khv',zpopthv(:,:,:,3))
1897      CALL send_xios_field('whv',zpopthv(:,:,:,4))
1898      CALL send_xios_field('ghv',zpopthv(:,:,:,5))
1899      IF (callclouds) THEN
1900         CALL send_xios_field('dtaui',zpoptti(:,:,:,1))
1901         CALL send_xios_field('taui',zpoptti(:,:,:,2))
1902         CALL send_xios_field('ki',zpoptti(:,:,:,3))
1903         CALL send_xios_field('wi',zpoptti(:,:,:,4))
1904         CALL send_xios_field('gi',zpoptti(:,:,:,5))
1905         CALL send_xios_field('dtauv',zpopttv(:,:,:,1))
1906         CALL send_xios_field('tauv',zpopttv(:,:,:,2))
1907         CALL send_xios_field('kv',zpopttv(:,:,:,3))
1908         CALL send_xios_field('wv',zpopttv(:,:,:,4))
1909         CALL send_xios_field('gv',zpopttv(:,:,:,5))
1910      ENDIF
1911
1912      !--------------------------------------------------------
1913      ! Microphysical tracers :
1914      !--------------------------------------------------------
1915      IF (callmufi) THEN
1916         ! Atmosphere (3D) :
1917         ! Moments M0 and M3 :
1918         CALL send_xios_field("mu_m0as",zq(:,:,micro_indx(1))*i2e)
1919         CALL send_xios_field("mu_m3as",zq(:,:,micro_indx(2))*i2e)
1920         CALL send_xios_field("mu_m0af",zq(:,:,micro_indx(3))*i2e)
1921         CALL send_xios_field("mu_m3af",zq(:,:,micro_indx(4))*i2e)
1922         IF (callclouds) THEN
1923            CALL send_xios_field("mu_m0n",zq(:,:,micro_indx(5))*i2e)
1924            CALL send_xios_field("mu_m3n",zq(:,:,micro_indx(6))*i2e)
1925            DO iq = 1, size(ices_indx)
1926               CALL send_xios_field(TRIM(nameOfTracer(ices_indx(iq))),zq(:,:,ices_indx(iq))*i2e)
1927            ENDDO
1928         ENDIF
1929
1930         ! Microphysical diagnostics :
1931         CALL send_xios_field("rc_sph",mmd_rc_sph(:,:))
1932         CALL send_xios_field("rc_fra",mmd_rc_fra(:,:))
1933         CALL send_xios_field("vsed_aers",mmd_aer_s_w(:,:))
1934         CALL send_xios_field("vsed_aerf",mmd_aer_f_w(:,:))
1935         CALL send_xios_field("flux_aers",mmd_aer_s_flux(:,:))
1936         CALL send_xios_field("flux_aerf",mmd_aer_f_flux(:,:))
1937         IF (callclouds) THEN
1938            CALL send_xios_field("rc_cld",mmd_rc_cld(:,:))
1939            CALL send_xios_field("vsed_ccn",mmd_ccn_w(:,:))
1940            CALL send_xios_field("flux_ccn",mmd_ccn_flux(:,:))
1941            DO iq = 1, size(ices_indx)
1942               CALL send_xios_field('flux_i'//TRIM(nameOfTracer(gazs_indx(iq))),mmd_ice_fluxes(:,:,iq))
1943               CALL send_xios_field(TRIM(nameOfTracer(gazs_indx(iq)))//'_sat',mmd_gazs_sat(:,:,iq))
1944            ENDDO
1945         ENDIF
1946         
1947         ! Surface (2D) :
1948         CALL send_xios_field("aer_prec",mmd_aer_prec(:))
1949         IF (callclouds) THEN
1950            CALL send_xios_field("ccn_prec",mmd_ccn_prec(:))
1951            DO iq = 1, size(ices_indx)
1952               CALL send_xios_field('i'//TRIM(nameOfTracer(gazs_indx(iq)))//'_prec',mmd_ice_prec(:,iq))
1953            ENDDO
1954         ENDIF
1955      ENDIF ! of 'if callmufi'
1956
1957      !--------------------------------------------------------
1958      ! Chemical tracers :
1959      !--------------------------------------------------------
1960      IF (callchim) THEN
1961         ! Surface (2D) :
1962         CALL send_xios_field("evapCH4",dycevapCH4(:)) ! Pseudo-evaporation flux (mol/mol/s)
1963         CALL send_xios_field("tankCH4",tankCH4(:))    ! CH4 tank at the surface (m)
1964
1965         ! Atmosphere (3D) :
1966         ! Chemical species :
1967         DO iq = 1, nkim
1968            ! If no cloud : gzs_indx uninitialized
1969            CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
1970         ENDDO
1971
1972         ! Condensation tendencies from microphysics (mol/mol/s) :
1973         IF (callclouds) THEN
1974            DO iq = 1, size(ices_indx)
1975               CALL send_xios_field('dmuficond_'//trim(nameOfTracer(gazs_indx(iq))),dmuficond(:,:,iq)/rat_mmol(gazs_indx(iq))) ! kg/kg/s -> mol/mol/s
1976            ENDDO
1977         ENDIF
1978         
1979         ! Condensation tendencies (mol/mol/s) :
1980         CALL send_xios_field("dqcond_CH4",dyccond(:,:,7+nmicro))
1981         CALL send_xios_field("dqcond_C2H2",dyccond(:,:,10+nmicro))
1982         CALL send_xios_field("dqcond_C2H4",dyccond(:,:,12+nmicro))
1983         CALL send_xios_field("dqcond_C2H6",dyccond(:,:,14+nmicro))
1984         CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro))
1985         CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro))
1986         CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro))
1987         CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro))
1988         CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro))
1989         CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro))
1990         CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro))
1991         CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro))
1992         CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro))
1993         CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro))
1994         CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro))
1995         CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro))
1996         CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro))
1997     
1998         ! Upper atmosphere chemistry variables (3D) :
1999         DO iq = 1, nkim
2000            CALL send_xios_field(trim(cnames(iq))//"_up",ykim_up(iq,:,:)) ! mol/mol
2001         ENDDO
2002         
2003         ! Total atmosphere chemistry variables (3D) :
2004         ! Append fields in ykim_tot for output on the total vertical grid (0->1300km)
2005         DO iq = 1, nkim   
2006            ! GCM levels
2007            DO l = 1, nlayer
2008               ykim_tot(iq,:,l) = zq(:,l,iq+nmicro)/rat_mmol(iq+nmicro)
2009            ENDDO
2010            ! Upper levels
2011            DO l = 1, nlaykim_up
2012               ykim_tot(iq,:,nlayer+l) = ykim_up(iq,:,l)
2013            ENDDO
2014            CALL send_xios_field(trim(cnames(iq))//"_tot",ykim_tot(iq,:,:)) ! mol/mol
2015         ENDDO
2016      ENDIF ! of 'if callchim'
2017
2018
2019      if (lastcall.and.is_omp_master) then
2020        write(*,*) "physiq: call xios_context_finalize"
2021        call xios_context_finalize
2022      endif
2023#endif
2024#else
2025      !MESOSCALE outputs     
2026      comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
2027      comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
2028      comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
2029      comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
2030      comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
2031      comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
2032      comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
2033      comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
2034      sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid)
2035      comm_zqfi_omp(1:ngrid,1:nlayer,1:nq) = zq(1:ngrid,1:nlayer,1:nq)
2036      comm_zdtlc(1:ngrid,1:nlayer) = zdtlc(1:ngrid,1:nlayer)
2037#endif     
2038
2039      icount=icount+1
2040
2041    end subroutine physiq
2042   
2043    end module physiq_mod
Note: See TracBrowser for help on using the repository browser.