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

Last change on this file since 1838 was 1838, checked in by mlefevre, 7 years ago

MESOSCALE GENERIC. Some declaration of allocatable remained in physiq_mod. Now done in phys_state_var_mod

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