source: trunk/LMDZ.MARS/libf/phymars/physiq.F @ 1307

Last change on this file since 1307 was 1292, checked in by aslmd, 11 years ago

LMDZ.MARS. MESOSCALE only. modif to output tau_ice when clouds not radiatively active.

File size: 98.1 KB
Line 
1      SUBROUTINE physiq(
2     $            ngrid,nlayer,nq
3     $            ,firstcall,lastcall
4     $            ,pday,ptime,ptimestep
5     $            ,pplev,pplay,pphi
6     $            ,pu,pv,pt,pq
7     $            ,pw
8     $            ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
9
10      use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2,
11     &                      igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
12     &                      igcm_ccn_mass, igcm_ccn_number,
13     &                      igcm_dust_mass, igcm_dust_number, igcm_h2o2,
14     &                      nuice_ref, rho_ice, rho_dust, ref_r0
15      use comsoil_h, only: inertiedat, ! soil thermal inertia
16     &                     tsoil, nsoilmx ! number of subsurface layers
17      use comgeomfi_h, only: long, lati, area,
18     &                       sinlon, coslon, sinlat, coslat
19      use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
20     &                     zthe, z0, albedo_h2o_ice,
21     &                     frost_albedo_threshold,
22     &                     tsurf, co2ice, emis,
23     &                     capcal, fluxgrd, qsurf
24      use comsaison_h, only: dist_sol, declin, mu0, fract
25      use slope_mod, only: theta_sl, psi_sl
26      use conc_mod, only: rnew, cpnew, mmean
27      use control_mod, only: iphysiq, day_step, ecritstart
28      use dimradmars_mod, only: tauscaling, aerosol,
29     &                          dtrad, fluxrad_sky, fluxrad, albedo,
30     &                          naerkind
31      use turb_mod, only: q2, wstar, ustar, sensibFlux,
32     &                    zmax_th, hfmax_th, turb_resolved
33      use planete_h, only: aphelie, periheli, year_day, peri_day,
34     &                     obliquit
35      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad, daysec
36      use param_v4_h, only: nreact,n_avog,
37     &                      fill_data_thermos, allocate_param_thermos
38      use iono_h, only: allocate_param_iono
39#ifdef MESOSCALE
40      use comsoil_h, only: mlayer,layer
41      use surfdat_h, only: z0_default
42      use comm_wrf
43#else
44      use phyredem, only: physdem0, physdem1
45      use eofdump_mod, only: eofdump
46#endif
47
48      IMPLICIT NONE
49c=======================================================================
50c
51c   subject:
52c   --------
53c
54c   Organisation of the physical parametrisations of the LMD
55c   martian atmospheric general circulation model.
56c
57c   The GCM can be run without or with tracer transport
58c   depending on the value of Logical "tracer" in file  "callphys.def"
59c   Tracers may be water vapor, ice OR chemical species OR dust particles
60c
61c   SEE comments in initracer.F about numbering of tracer species...
62c
63c   It includes:
64c
65c      1. Initialization:
66c      1.1 First call initializations
67c      1.2 Initialization for every call to physiq
68c      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
69c      2. Compute radiative transfer tendencies
70c         (longwave and shortwave) for CO2 and aerosols.
71c      3. Gravity wave and subgrid scale topography drag :
72c      4. Vertical diffusion (turbulent mixing):
73c      5. Convective adjustment
74c      6. Condensation and sublimation of carbon dioxide.
75c      7.  TRACERS :
76c       7a. water and water ice
77c       7b. call for photochemistry when tracers are chemical species
78c       7c. other scheme for tracer (dust) transport (lifting, sedimentation)
79c       7d. updates (CO2 pressure variations, surface budget)
80c      8. Contribution to tendencies due to thermosphere
81c      9. Surface and sub-surface temperature calculations
82c     10. Write outputs :
83c           - "startfi", "histfi" (if it's time)
84c           - Saving statistics (if "callstats = .true.")
85c           - Dumping eof (if "calleofdump = .true.")
86c           - Output any needed variables in "diagfi"
87c     11. Diagnostic: mass conservation of tracers
88c
89c   author:
90c   -------
91c           Frederic Hourdin    15/10/93
92c           Francois Forget             1994
93c           Christophe Hourdin  02/1997
94c           Subroutine completly rewritten by F.Forget (01/2000)
95c           Introduction of the photochemical module: S. Lebonnois (11/2002)
96c           Introduction of the thermosphere module: M. Angelats i Coll (2002)
97c           Water ice clouds: Franck Montmessin (update 06/2003)
98c           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
99c             Nb: See callradite.F for more information.
100c           Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
101c           jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
102c           
103c   arguments:
104c   ----------
105c
106c   input:
107c   ------
108c    ecri                  period (in dynamical timestep) to write output
109c    ngrid                 Size of the horizontal grid.
110c                          All internal loops are performed on that grid.
111c    nlayer                Number of vertical layers.
112c    nq                    Number of advected fields
113c    firstcall             True at the first call
114c    lastcall              True at the last call
115c    pday                  Number of days counted from the North. Spring
116c                          equinoxe.
117c    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
118c    ptimestep             timestep (s)
119c    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
120c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
121c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
122c    pu(ngrid,nlayer)      u component of the wind (ms-1)
123c    pv(ngrid,nlayer)      v component of the wind (ms-1)
124c    pt(ngrid,nlayer)      Temperature (K)
125c    pq(ngrid,nlayer,nq)   Advected fields
126c    pudyn(ngrid,nlayer)    |
127c    pvdyn(ngrid,nlayer)    | Dynamical temporal derivative for the
128c    ptdyn(ngrid,nlayer)    | corresponding variables
129c    pqdyn(ngrid,nlayer,nq) |
130c    pw(ngrid,?)           vertical velocity
131c
132c   output:
133c   -------
134c
135c    pdu(ngrid,nlayer)       |
136c    pdv(ngrid,nlayer)       |  Temporal derivative of the corresponding
137c    pdt(ngrid,nlayer)       |  variables due to physical processes.
138c    pdq(ngrid,nlayer,nq)    |
139c    pdpsrf(ngrid)           |
140c    tracerdyn                 call tracer in dynamical part of GCM ?
141
142c
143c=======================================================================
144c
145c    0.  Declarations :
146c    ------------------
147
148#include "dimensions.h"
149#include "callkeys.h"
150#include "comg1d.h"
151#include "nlteparams.h"
152#include "comvert.h"
153#include "chimiedata.h"
154#include "netcdf.inc"
155
156c Arguments :
157c -----------
158
159c   inputs:
160c   -------
161      INTEGER ngrid,nlayer,nq
162      REAL ptimestep
163      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
164      REAL pphi(ngrid,nlayer)
165      REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
166      REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
167      REAL pw(ngrid,nlayer) !Mars pvervel transmit par dyn3d
168      REAL zh(ngrid,nlayer)      ! potential temperature (K)
169      LOGICAL firstcall,lastcall
170
171      REAL pday
172      REAL ptime
173      logical tracerdyn
174
175c   outputs:
176c   --------
177c     physical tendencies
178      REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)
179      REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
180      REAL pdpsrf(ngrid) ! surface pressure tendency
181
182
183c Local saved variables:
184c ----------------------
185      INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
186      INTEGER,SAVE :: icount     ! counter of calls to physiq during the run.
187     
188c     Variables used by the water ice microphysical scheme:
189      REAL rice(ngrid,nlayer)    ! Water ice geometric mean radius (m)
190      REAL nuice(ngrid,nlayer)   ! Estimated effective variance
191                                     !   of the size distribution
192      real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
193      real rhocloud(ngrid,nlayer)  ! Cloud density (kg.m-3)
194      REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
195      REAL surfice(ngrid,nlayer)  !  ice surface area (m2/m3, if photochemistry)
196      REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
197                                      ! thermal inertia (J.s-1/2.m-2.K-1)
198                                      ! (used only when tifeedback=.true.)
199
200c     Variables used by the slope model
201      REAL sl_ls, sl_lct, sl_lat
202      REAL sl_tau, sl_alb, sl_the, sl_psi
203      REAL sl_fl0, sl_flu
204      REAL sl_ra, sl_di0
205      REAL sky
206
207      REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant
208
209c Local variables :
210c -----------------
211
212      REAL CBRT
213      EXTERNAL CBRT
214
215!      CHARACTER*80 fichier
216      INTEGER l,ig,ierr,igout,iq,tapphys
217
218      REAL fluxsurf_lw(ngrid)      !incident LW (IR) surface flux (W.m-2)
219      REAL fluxsurf_sw(ngrid,2)    !incident SW (solar) surface flux (W.m-2)
220      REAL fluxtop_lw(ngrid)       !Outgoing LW (IR) flux to space (W.m-2)
221      REAL fluxtop_sw(ngrid,2)     !Outgoing SW (solar) flux to space (W.m-2)
222      REAL tauref(ngrid)           ! Reference column optical depth at odpref
223      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
224      REAL tau(ngrid,naerkind)     ! Column dust optical depth at each point
225                                   ! AS: TBD: this one should be in a module !
226      REAL zls                       !  solar longitude (rad)
227      REAL zday                      ! date (time since Ls=0, in martian days)
228      REAL zzlay(ngrid,nlayer)     ! altitude at the middle of the layers
229      REAL zzlev(ngrid,nlayer+1)   ! altitude at layer boundaries
230!      REAL latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
231
232c     Tendancies due to various processes:
233      REAL dqsurf(ngrid,nq)
234      REAL zdtlw(ngrid,nlayer)     ! (K/s)
235      REAL zdtsw(ngrid,nlayer)     ! (K/s)
236!      REAL cldtlw(ngrid,nlayer)     ! (K/s) LW heating rate for clear area
237!      REAL cldtsw(ngrid,nlayer)     ! (K/s) SW heating rate for clear area
238      REAL zdtnirco2(ngrid,nlayer) ! (K/s)
239      REAL zdtnlte(ngrid,nlayer)   ! (K/s)
240      REAL zdtsurf(ngrid)            ! (K/s)
241      REAL zdtcloud(ngrid,nlayer)
242      REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer)  ! (m.s-2)
243      REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid)         ! (K/s)
244      REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer)  ! (m.s-2)
245      REAL zdhadj(ngrid,nlayer)                           ! (K/s)
246      REAL zdtgw(ngrid,nlayer)                            ! (K/s)
247      REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer)    ! (m.s-2)
248      REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid)
249      REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
250
251      REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
252      REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
253      REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
254      REAL zdqadj(ngrid,nlayer,nq)
255      REAL zdqc(ngrid,nlayer,nq)
256      REAL zdqcloud(ngrid,nlayer,nq)
257      REAL zdqscloud(ngrid,nq)
258      REAL zdqchim(ngrid,nlayer,nq)
259      REAL zdqschim(ngrid,nq)
260
261      REAL zdteuv(ngrid,nlayer)    ! (K/s)
262      REAL zdtconduc(ngrid,nlayer) ! (K/s)
263      REAL zdumolvis(ngrid,nlayer)
264      REAL zdvmolvis(ngrid,nlayer)
265      real zdqmoldiff(ngrid,nlayer,nq)
266
267c     Local variable for local intermediate calcul:
268      REAL zflubid(ngrid)
269      REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
270      REAL zdum1(ngrid,nlayer)
271      REAL zdum2(ngrid,nlayer)
272      REAL ztim1,ztim2,ztim3, z1,z2
273      REAL ztime_fin
274      REAL zdh(ngrid,nlayer)
275      INTEGER length
276      PARAMETER (length=100)
277
278c local variables only used for diagnostic (output in file "diagfi" or "stats")
279c -----------------------------------------------------------------------------
280      REAL ps(ngrid), zt(ngrid,nlayer)
281      REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
282      REAL zq(ngrid,nlayer,nq)
283      REAL fluxtop_sw_tot(ngrid), fluxsurf_sw_tot(ngrid)
284      character*2 str2
285!      character*5 str5
286      real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
287      real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
288      integer igmin, lmin
289      logical tdiag
290
291      real co2col(ngrid)        ! CO2 column
292      ! pplev and pplay are dynamical inputs and must not be modified in the physics.
293      ! instead, use zplay and zplev :
294      REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
295!      REAL zstress(ngrid),cd
296      real tmean, zlocal(nlayer)
297      real rho(ngrid,nlayer)  ! density
298      real vmr(ngrid,nlayer)  ! volume mixing ratio
299      real rhopart(ngrid,nlayer) ! number density of a given species
300      real colden(ngrid,nq)     ! vertical column of tracers
301      REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
302      REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
303      REAL Nccntot(ngrid)       ! Total number of ccn (nbr/m2)
304      REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
305      REAL rave(ngrid)          ! Mean water ice effective radius (m)
306      REAL opTES(ngrid,nlayer)  ! abs optical depth at 825 cm-1
307      REAL tauTES(ngrid)        ! column optical depth at 825 cm-1
308      REAL Qabsice                ! Water ice absorption coefficient
309      REAL taucloudtes(ngrid)  ! Cloud opacity at infrared
310                               !   reference wavelength using
311                               !   Qabs instead of Qext
312                               !   (direct comparison with TES)
313                               
314      REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
315      REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
316      REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
317      REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
318      REAL nccn(ngrid,nlayer)  ! true n ccn (kg/kg)
319      REAL qccn(ngrid,nlayer)  ! true q ccn (kg/kg)
320
321c Test 1d/3d scavenging
322      real h2otot(ngrid)
323      REAL satu(ngrid,nlayer)  ! satu ratio for output
324      REAL zqsat(ngrid,nlayer) ! saturation
325
326      REAL,SAVE :: time_phys
327
328! Added for new NLTE scheme
329
330      real co2vmr_gcm(ngrid,nlayer)
331      real n2vmr_gcm(ngrid,nlayer)
332      real ovmr_gcm(ngrid,nlayer)
333      real covmr_gcm(ngrid,nlayer)
334      integer ierr_nlte
335      real*8  varerr
336
337c Variables for PBL
338      REAL zz1(ngrid)
339      REAL lmax_th_out(ngrid)
340      REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
341      REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
342      INTEGER lmax_th(ngrid),dimout,n_out,n
343      CHARACTER(50) zstring
344      REAL dtke_th(ngrid,nlayer+1)
345      REAL zcdv(ngrid), zcdh(ngrid)
346      REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
347      REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
348      REAL u_out1(ngrid)
349      REAL T_out1(ngrid)
350      REAL, ALLOCATABLE, DIMENSION(:) :: z_out     ! height of interpolation between z0 and z1 [meters]
351      REAL tstar(ngrid)  ! friction velocity and friction potential temp
352      REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
353!      REAL zu2(ngrid)
354
355c=======================================================================
356
357c 1. Initialisation:
358c -----------------
359
360c  1.1   Initialisation only at first call
361c  ---------------------------------------
362      IF (firstcall) THEN
363
364c        variables set to 0
365c        ~~~~~~~~~~~~~~~~~~
366         aerosol(:,:,:)=0
367         dtrad(:,:)=0
368
369#ifndef MESOSCALE
370         fluxrad(:)=0
371         wstar(:)=0.
372#endif
373
374c        read startfi
375c        ~~~~~~~~~~~~
376#ifndef MESOSCALE
377! GCM. Read netcdf initial physical parameters.
378         CALL phyetat0 ("startfi.nc",0,0,
379     &         nsoilmx,ngrid,nlayer,nq,
380     &         day_ini,time_phys,
381     &         tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)
382
383         if (pday.ne.day_ini) then
384           write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
385     &                "physics and dynamics"
386           write(*,*) "dynamics day: ",pday
387           write(*,*) "physics day:  ",day_ini
388           stop
389         endif
390
391         write (*,*) 'In physiq day_ini =', day_ini
392
393#else
394! MESOSCALE. Supposedly everything is already set in modules.
395! So we just check. And we calculate day_ini + two param in dyn3d/control_mod.
396      print*,"check: rad,cpp,g,r,rcp,daysec"
397      print*,rad,cpp,g,r,rcp,daysec
398      PRINT*,'check: tsurf ',tsurf(1),tsurf(ngrid)
399      PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngrid,nsoilmx)
400      PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
401      PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
402      PRINT*,'check: tracernames ', noms
403      PRINT*,'check: emis ',emis(1),emis(ngrid)
404      PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
405      PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq)
406      PRINT*,'check: co2 ',co2ice(1),co2ice(ngrid)
407      day_step=daysec/ptimestep
408      PRINT*,'Call to LMD physics:',day_step,' per Martian day'
409      iphysiq=ptimestep
410      day_ini = pday
411#endif
412
413c        initialize tracers
414c        ~~~~~~~~~~~~~~~~~~
415         tracerdyn=tracer
416         IF (tracer) THEN
417            CALL initracer(ngrid,nq,qsurf)
418         ENDIF  ! end tracer
419
420c        Initialize albedo and orbital calculation
421c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
422         CALL surfini(ngrid,co2ice,qsurf,albedo)
423         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
424
425c        initialize soil
426c        ~~~~~~~~~~~~~~~
427         IF (callsoil) THEN
428c           Thermal inertia feedback:
429            IF (tifeedback) THEN
430                CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
431                CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
432     s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
433            ELSE
434                CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
435     s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
436            ENDIF ! of IF (tifeedback)
437         ELSE
438            PRINT*,
439     &     'PHYSIQ WARNING! Thermal conduction in the soil turned off'
440            DO ig=1,ngrid
441               capcal(ig)=1.e5
442               fluxgrd(ig)=0.
443            ENDDO
444         ENDIF
445         icount=1
446
447#ifndef MESOSCALE
448c        Initialize thermospheric parameters
449c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
450
451         if (callthermos) then
452            call fill_data_thermos
453            call allocate_param_thermos(nlayer)
454            call allocate_param_iono(nlayer,nreact)
455            if(solvarmod.eq.0) call param_read
456            if(solvarmod.eq.1) call param_read_e107
457         endif
458#endif
459c        Initialize R and Cp as constant
460
461         if (.not.callthermos .and. .not.photochem) then
462                 do l=1,nlayer
463                  do ig=1,ngrid
464                   rnew(ig,l)=r
465                   cpnew(ig,l)=cpp
466                   mmean(ig,l)=mugaz
467                   enddo
468                  enddo 
469         endif         
470
471         if(callnlte.and.nltemodel.eq.2) call nlte_setup
472         if(callnirco2.and.nircorr.eq.1) call NIR_leedat
473         if(thermochem) call chemthermos_readini
474
475        IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
476          write(*,*)"physiq: water_param Surface water ice albedo:",
477     .                  albedo_h2o_ice
478        ENDIF
479
480#ifndef MESOSCALE
481         if (callslope) call getslopes(ngrid,phisfi)
482
483         if (ngrid.ne.1) then ! no need to create a restart file in 1d
484         call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq,
485     .                 ptimestep,pday,time_phys,area,
486     .                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
487         endif
488#endif
489                   
490      ENDIF        !  (end of "if firstcall")
491
492
493c ---------------------------------------------------
494c 1.2   Initializations done at every physical timestep:
495c ---------------------------------------------------
496c
497
498c     Initialize various variables
499c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
500      pdv(:,:)=0
501      pdu(:,:)=0
502      pdt(:,:)=0
503      pdq(:,:,:)=0
504      pdpsrf(:)=0
505      zflubid(:)=0
506      zdtsurf(:)=0
507      dqsurf(:,:)=0
508      igout=ngrid/2+1
509
510
511      zday=pday+ptime ! compute time, in sols (and fraction thereof)
512
513c     Compute Solar Longitude (Ls) :
514c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
515      if (season) then
516         call solarlong(zday,zls)
517      else
518         call solarlong(float(day_ini),zls)
519      end if
520
521c     Initialize pressure levels
522c     ~~~~~~~~~~~~~~~~~~
523      zplev(:,:) = pplev(:,:)
524      zplay(:,:) = pplay(:,:)
525      ps(:) = pplev(:,1)
526
527
528c     Compute geopotential at interlayers
529c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
530c     ponderation des altitudes au niveau des couches en dp/p
531
532      DO l=1,nlayer
533         DO ig=1,ngrid
534            zzlay(ig,l)=pphi(ig,l)/g
535         ENDDO
536      ENDDO
537      DO ig=1,ngrid
538         zzlev(ig,1)=0.
539         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
540      ENDDO
541      DO l=2,nlayer
542         DO ig=1,ngrid
543            z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
544            z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
545            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
546         ENDDO
547      ENDDO
548
549
550!     Potential temperature calculation not the same in physiq and dynamic
551
552c     Compute potential temperature
553c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
554      DO l=1,nlayer
555         DO ig=1,ngrid
556            zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
557            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
558         ENDDO
559      ENDDO
560
561#ifndef MESOSCALE
562c-----------------------------------------------------------------------
563c    1.2.5 Compute mean mass, cp, and R
564c    --------------------------------
565
566      if(photochem.or.callthermos) then
567         call concentrations(ngrid,nlayer,nq,
568     &                       zplay,pt,pdt,pq,pdq,ptimestep)
569      endif
570#endif
571c-----------------------------------------------------------------------
572c    2. Compute radiative tendencies :
573c------------------------------------
574
575
576      IF (callrad) THEN
577         IF( MOD(icount-1,iradia).EQ.0) THEN
578
579c          Local Solar zenith angle
580c          ~~~~~~~~~~~~~~~~~~~~~~~~
581           CALL orbite(zls,dist_sol,declin)
582
583           IF(diurnal) THEN
584               ztim1=SIN(declin)
585               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
586               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
587
588               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
589     s         ztim1,ztim2,ztim3, mu0,fract)
590
591           ELSE
592               CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
593           ENDIF
594
595c          NLTE cooling from CO2 emission
596c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
597           IF(callnlte) then
598              if(nltemodel.eq.0.or.nltemodel.eq.1) then
599                 CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
600              else if(nltemodel.eq.2) then
601                co2vmr_gcm(1:ngrid,1:nlayer)=
602     &                      pq(1:ngrid,1:nlayer,igcm_co2)*
603     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
604                n2vmr_gcm(1:ngrid,1:nlayer)=
605     &                      pq(1:ngrid,1:nlayer,igcm_n2)*
606     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
607                covmr_gcm(1:ngrid,1:nlayer)=
608     &                      pq(1:ngrid,1:nlayer,igcm_co)*
609     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
610                ovmr_gcm(1:ngrid,1:nlayer)=
611     &                      pq(1:ngrid,1:nlayer,igcm_o)*
612     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
613                 
614                 CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
615     $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
616     $                ovmr_gcm,  zdtnlte,ierr_nlte,varerr )
617                 if(ierr_nlte.gt.0) then
618                    write(*,*)
619     $                'WARNING: nlte_tcool output with error message',
620     $                'ierr_nlte=',ierr_nlte,'varerr=',varerr
621                    write(*,*)'I will continue anyway'
622                 endif
623
624                 zdtnlte(1:ngrid,1:nlayer)=
625     &                             zdtnlte(1:ngrid,1:nlayer)/86400.
626              endif
627           else
628              zdtnlte(:,:)=0.
629           endif
630
631c          Find number of layers for LTE radiation calculations
632           IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
633     &          CALL nlthermeq(ngrid,nlayer,zplev,zplay)
634
635c          Note: Dustopacity.F has been transferred to callradite.F
636         
637c          Call main radiative transfer scheme
638c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
639c          Transfer through CO2 (except NIR CO2 absorption)
640c            and aerosols (dust and water ice)
641
642c          Radiative transfer
643c          ------------------
644           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
645     $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
646     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
647     $     tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
648     $     nuice,co2ice)
649
650c          Outputs for basic check (middle of domain)
651c          ------------------------------------------
652           write(*,'("Ls =",f11.6," check lat =",f10.6,
653     &               " lon =",f11.6)')
654     &           zls*180./pi,lati(igout)*180/pi,long(igout)*180/pi
655           write(*,'(" tauref(",f4.0," Pa) =",f9.6,
656     &             " tau(",f4.0," Pa) =",f9.6)')
657     &            odpref,tauref(igout),
658     &            odpref,tau(igout,1)*odpref/zplev(igout,1)
659c          ---------------------------------------------------------
660c          Call slope parameterization for direct and scattered flux
661c          ---------------------------------------------------------
662           IF(callslope) THEN
663            print *, 'Slope scheme is on and computing...'
664            DO ig=1,ngrid 
665              sl_the = theta_sl(ig)
666              IF (sl_the .ne. 0.) THEN
667                ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
668                DO l=1,2
669                 sl_lct = ptime*24. + 180.*long(ig)/pi/15.
670                 sl_ra  = pi*(1.0-sl_lct/12.)
671                 sl_lat = 180.*lati(ig)/pi
672                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
673                 sl_alb = albedo(ig,l)
674                 sl_psi = psi_sl(ig)
675                 sl_fl0 = fluxsurf_sw(ig,l)
676                 sl_di0 = 0.
677                 if (mu0(ig) .gt. 0.) then
678                  sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
679                  sl_di0 = sl_di0*1370./dist_sol/dist_sol
680                  sl_di0 = sl_di0/ztim1
681                  sl_di0 = fluxsurf_sw(ig,l)*sl_di0
682                 endif
683                 ! you never know (roundup concern...)
684                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
685                 !!!!!!!!!!!!!!!!!!!!!!!!!!
686                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
687     &                             sl_tau, sl_alb, sl_the, sl_psi,
688     &                             sl_di0, sl_fl0, sl_flu )
689                 !!!!!!!!!!!!!!!!!!!!!!!!!!
690                 fluxsurf_sw(ig,l) = sl_flu
691                ENDDO
692              !!! compute correction on IR flux as well
693              sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
694              fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
695              ENDIF
696            ENDDO
697           ENDIF
698
699c          CO2 near infrared absorption
700c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
701           zdtnirco2(:,:)=0
702           if (callnirco2) then
703              call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
704     .                       mu0,fract,declin, zdtnirco2)
705           endif
706
707c          Radiative flux from the sky absorbed by the surface (W.m-2)
708c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
709           DO ig=1,ngrid
710               fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
711     $         +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
712     $         +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
713           ENDDO
714
715
716c          Net atmospheric radiative heating rate (K.s-1)
717c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
718           IF(callnlte) THEN
719              CALL blendrad(ngrid, nlayer, zplay,
720     &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
721           ELSE
722              DO l=1,nlayer
723                 DO ig=1,ngrid
724                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
725     &                          +zdtnirco2(ig,l)
726                  ENDDO
727              ENDDO
728           ENDIF
729
730        ENDIF ! of if(mod(icount-1,iradia).eq.0)
731
732c    Transformation of the radiative tendencies:
733c    -------------------------------------------
734
735c          Net radiative surface flux (W.m-2)
736c          ~~~~~~~~~~~~~~~~~~~~~~~~~~
737c
738           DO ig=1,ngrid
739               zplanck(ig)=tsurf(ig)*tsurf(ig)
740               zplanck(ig)=emis(ig)*
741     $         stephan*zplanck(ig)*zplanck(ig)
742               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
743               IF(callslope) THEN
744                 sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
745                 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
746               ENDIF
747           ENDDO
748
749         DO l=1,nlayer
750            DO ig=1,ngrid
751               pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
752            ENDDO
753         ENDDO
754
755      ENDIF ! of IF (callrad)
756
757c-----------------------------------------------------------------------
758c    3. Gravity wave and subgrid scale topography drag :
759c    -------------------------------------------------
760
761
762      IF(calllott)THEN
763
764        CALL calldrag_noro(ngrid,nlayer,ptimestep,
765     &                 zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
766 
767        DO l=1,nlayer
768          DO ig=1,ngrid
769            pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
770            pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
771            pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
772          ENDDO
773        ENDDO
774      ENDIF
775
776c-----------------------------------------------------------------------
777c    4. Vertical diffusion (turbulent mixing):
778c    -----------------------------------------
779
780      IF (calldifv) THEN
781
782         DO ig=1,ngrid
783            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
784         ENDDO
785
786         zdum1(:,:)=0
787         zdum2(:,:)=0
788         do l=1,nlayer
789            do ig=1,ngrid
790               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
791            enddo
792         enddo
793
794c ----------------------
795c Treatment of a special case : using new surface layer (Richardson based)
796c without using the thermals in gcm and mesoscale can yield problems in
797c weakly unstable situations when winds are near to 0. For those cases, we add
798c a unit subgrid gustiness. Remember that thermals should be used we using the
799c Richardson based surface layer model.
800        IF ( .not.calltherm
801     .       .and. callrichsl
802     .       .and. .not.turb_resolved) THEN
803          DO ig=1, ngrid
804             IF (zh(ig,1) .lt. tsurf(ig)) THEN
805               wstar(ig)=1.
806               hfmax_th(ig)=0.2
807             ELSE
808               wstar(ig)=0.
809               hfmax_th(ig)=0.
810             ENDIF     
811          ENDDO
812        ENDIF
813c ----------------------
814
815         IF (tke_heat_flux .ne. 0.) THEN
816             zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)*
817     &                 (-alog(zplay(:,1)/zplev(:,1)))
818             pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
819         ENDIF
820
821c        Calling vdif (Martian version WITH CO2 condensation)
822         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
823     $        ptimestep,capcal,lwrite,
824     $        zplay,zplev,zzlay,zzlev,z0,
825     $        pu,pv,zh,pq,tsurf,emis,qsurf,
826     $        zdum1,zdum2,zdh,pdq,zflubid,
827     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
828     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux)
829
830
831          DO ig=1,ngrid
832             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
833          ENDDO
834
835         IF (.not.turb_resolved) THEN
836          DO l=1,nlayer
837            DO ig=1,ngrid
838               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
839               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
840               pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
841
842               zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
843            ENDDO
844          ENDDO
845
846          if (tracer) then
847           DO iq=1, nq
848            DO l=1,nlayer
849              DO ig=1,ngrid
850                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
851              ENDDO
852            ENDDO
853           ENDDO
854           DO iq=1, nq
855              DO ig=1,ngrid
856                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
857              ENDDO
858           ENDDO
859          end if ! of if (tracer)
860         ELSE
861           write (*,*) '******************************************'
862           write (*,*) '** LES mode: the difv part is only used to'
863           write (*,*) '**  - provide HFX and UST to the dynamics'
864           write (*,*) '**  - update TSURF'
865           write (*,*) '******************************************'
866           !! Specific treatment for lifting in turbulent-resolving mode (AC)
867           IF (lifting .and. doubleq) THEN
868             !! lifted dust is injected in the first layer.
869             !! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
870             !! => lifted dust is not incremented before the sedimentation step.
871             zdqdif(1:ngrid,1,1:nq)=0.
872             zdqdif(1:ngrid,1,igcm_dust_number) =
873     .                  -zdqsdif(1:ngrid,igcm_dust_number)
874             zdqdif(1:ngrid,1,igcm_dust_mass) =
875     .                  -zdqsdif(1:ngrid,igcm_dust_mass)
876             zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
877             DO iq=1, nq
878              IF ((iq .ne. igcm_dust_mass)
879     &        .and. (iq .ne. igcm_dust_number)) THEN
880                zdqsdif(:,iq)=0.
881              ENDIF
882             ENDDO
883           ELSE
884             zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
885             zdqsdif(1:ngrid,1:nq) = 0.
886           ENDIF
887         ENDIF
888      ELSE   
889         DO ig=1,ngrid
890            zdtsurf(ig)=zdtsurf(ig)+
891     s        (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
892         ENDDO
893         IF (turb_resolved) THEN
894            write(*,*) 'Turbulent-resolving mode !'
895            write(*,*) 'Please set calldifv to T in callphys.def'
896            STOP
897         ENDIF
898      ENDIF ! of IF (calldifv)
899
900c-----------------------------------------------------------------------
901c   5. Thermals :
902c   -----------------------------
903
904      if(calltherm .and. .not.turb_resolved) then
905 
906        call calltherm_interface(ngrid,nlayer,nq,
907     $ tracer,igcm_co2,
908     $ zzlev,zzlay,
909     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
910     $ zplay,zplev,pphi,zpopsk,
911     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
912     $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
913     
914         DO l=1,nlayer
915           DO ig=1,ngrid
916              pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
917              pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
918              pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
919              q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
920           ENDDO
921        ENDDO
922 
923        DO ig=1,ngrid
924          q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
925        ENDDO     
926   
927        if (tracer) then
928        DO iq=1,nq
929         DO l=1,nlayer
930           DO ig=1,ngrid
931             pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
932           ENDDO
933         ENDDO
934        ENDDO
935        endif
936
937        lmax_th_out(:)=real(lmax_th(:))
938
939      else   !of if calltherm
940        lmax_th(:)=0
941        wstar(:)=0.
942        hfmax_th(:)=0.
943        lmax_th_out(:)=0.
944      end if
945
946c-----------------------------------------------------------------------
947c   5. Dry convective adjustment:
948c   -----------------------------
949
950      IF(calladj) THEN
951
952         DO l=1,nlayer
953            DO ig=1,ngrid
954               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
955            ENDDO
956         ENDDO
957         zduadj(:,:)=0
958         zdvadj(:,:)=0
959         zdhadj(:,:)=0
960
961         CALL convadj(ngrid,nlayer,nq,ptimestep,
962     $                zplay,zplev,zpopsk,lmax_th,
963     $                pu,pv,zh,pq,
964     $                pdu,pdv,zdh,pdq,
965     $                zduadj,zdvadj,zdhadj,
966     $                zdqadj)
967
968
969         DO l=1,nlayer
970            DO ig=1,ngrid
971               pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
972               pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
973               pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
974
975               zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
976            ENDDO
977         ENDDO
978
979         if(tracer) then
980           DO iq=1, nq
981            DO l=1,nlayer
982              DO ig=1,ngrid
983                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
984              ENDDO
985            ENDDO
986           ENDDO
987         end if
988      ENDIF ! of IF(calladj)
989
990
991
992c-----------------------------------------------------------------------
993c   6. Specific parameterizations for tracers
994c:   -----------------------------------------
995
996      if (tracer) then
997
998c   6a. Water and ice
999c     ---------------
1000
1001c        ---------------------------------------
1002c        Water ice condensation in the atmosphere
1003c        ----------------------------------------
1004         IF (water) THEN
1005
1006           call watercloud(ngrid,nlayer,ptimestep,
1007     &                zplev,zplay,pdpsrf,zzlay, pt,pdt,
1008     &                pq,pdq,zdqcloud,zdtcloud,
1009     &                nq,tau,tauscaling,rdust,rice,nuice,
1010     &                rsedcloud,rhocloud)
1011     
1012c Temperature variation due to latent heat release
1013           if (activice) then
1014               pdt(1:ngrid,1:nlayer) =
1015     &          pdt(1:ngrid,1:nlayer) +
1016     &          zdtcloud(1:ngrid,1:nlayer)
1017           endif
1018           
1019! increment water vapour and ice atmospheric tracers tendencies
1020           pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
1021     &       pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
1022     &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
1023           pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
1024     &       pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
1025     &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
1026
1027! increment dust and ccn masses and numbers
1028! We need to check that we have Nccn & Ndust > 0
1029! This is due to single precision rounding problems
1030           if (microphys) then
1031             pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
1032     &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
1033     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
1034             pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
1035     &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
1036     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
1037             where (pq(:,:,igcm_ccn_mass) +
1038     &       ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
1039               pdq(:,:,igcm_ccn_mass) =
1040     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1041               pdq(:,:,igcm_ccn_number) =
1042     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1043             end where
1044             where (pq(:,:,igcm_ccn_number) +
1045     &       ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
1046               pdq(:,:,igcm_ccn_mass) =
1047     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
1048               pdq(:,:,igcm_ccn_number) =
1049     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
1050             end where
1051           endif
1052
1053           if (scavenging) then
1054             pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
1055     &         pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
1056     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
1057             pdq(1:ngrid,1:nlayer,igcm_dust_number) =
1058     &         pdq(1:ngrid,1:nlayer,igcm_dust_number) +
1059     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
1060             where (pq(:,:,igcm_dust_mass) +
1061     &       ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
1062               pdq(:,:,igcm_dust_mass) =
1063     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1064               pdq(:,:,igcm_dust_number) =
1065     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1066             end where
1067             where (pq(:,:,igcm_dust_number) +
1068     &       ptimestep*pdq(:,:,igcm_dust_number) < 0.)
1069               pdq(:,:,igcm_dust_mass) =
1070     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
1071               pdq(:,:,igcm_dust_number) =
1072     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
1073             end where
1074           endif ! of if scavenging
1075                     
1076
1077         END IF  ! of IF (water)
1078
1079c   6b. Aerosol particles
1080c     -------------------
1081
1082c        ----------
1083c        Dust devil :
1084c        ----------
1085         IF(callddevil) then
1086           call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
1087     &                zdqdev,zdqsdev)
1088 
1089           if (dustbin.ge.1) then
1090              do iq=1,nq
1091                 DO l=1,nlayer
1092                    DO ig=1,ngrid
1093                       pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
1094                    ENDDO
1095                 ENDDO
1096              enddo
1097              do iq=1,nq
1098                 DO ig=1,ngrid
1099                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
1100                 ENDDO
1101              enddo
1102           endif  ! of if (dustbin.ge.1)
1103
1104         END IF ! of IF (callddevil)
1105
1106c        -------------
1107c        Sedimentation :   acts also on water ice
1108c        -------------
1109         IF (sedimentation) THEN
1110           zdqsed(1:ngrid,1:nlayer,1:nq)=0
1111           zdqssed(1:ngrid,1:nq)=0
1112
1113           call callsedim(ngrid,nlayer, ptimestep,
1114     &                zplev,zzlev, zzlay, pt, pdt, rdust, rice,
1115     &                rsedcloud,rhocloud,
1116     &                pq, pdq, zdqsed, zdqssed,nq,
1117     &                tau,tauscaling)
1118
1119           DO iq=1, nq
1120             DO l=1,nlayer
1121               DO ig=1,ngrid
1122                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
1123               ENDDO
1124             ENDDO
1125           ENDDO
1126           DO iq=1, nq
1127             DO ig=1,ngrid
1128                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
1129             ENDDO
1130           ENDDO
1131         END IF   ! of IF (sedimentation)
1132       
1133c Add lifted dust to tendancies after sedimentation in the LES (AC)
1134      IF (turb_resolved) THEN
1135           DO iq=1, nq
1136            DO l=1,nlayer
1137              DO ig=1,ngrid
1138                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
1139              ENDDO
1140            ENDDO
1141           ENDDO
1142           DO iq=1, nq
1143              DO ig=1,ngrid
1144                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
1145              ENDDO
1146           ENDDO
1147      ENDIF
1148
1149c
1150c   6c. Chemical species
1151c     ------------------
1152
1153#ifndef MESOSCALE
1154c        --------------
1155c        photochemistry :
1156c        --------------
1157         IF (photochem .or. thermochem) then
1158
1159!           dust and ice surface area
1160            call surfacearea(ngrid, nlayer, naerkind,
1161     $                       ptimestep, zplay, zzlay,
1162     $                       pt, pq, pdq, nq,
1163     $                       rdust, rice, tau, tauscaling,
1164     $                       surfdust, surfice)
1165!           call photochemistry
1166            call calchim(ngrid,nlayer,nq,
1167     &                   ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
1168     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
1169     $                   zdqcloud,zdqscloud,tauref,co2ice,
1170     $                   pu,pdu,pv,pdv,surfdust,surfice)
1171
1172           ! increment values of tracers:
1173           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1174                      ! tracers is zero anyways
1175             DO l=1,nlayer
1176               DO ig=1,ngrid
1177                 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
1178               ENDDO
1179             ENDDO
1180           ENDDO ! of DO iq=1,nq
1181           
1182           ! add condensation tendency for H2O2
1183           if (igcm_h2o2.ne.0) then
1184             DO l=1,nlayer
1185               DO ig=1,ngrid
1186                 pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
1187     &                                +zdqcloud(ig,l,igcm_h2o2)
1188               ENDDO
1189             ENDDO
1190           endif
1191
1192           ! increment surface values of tracers:
1193           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
1194                      ! tracers is zero anyways
1195             DO ig=1,ngrid
1196               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
1197             ENDDO
1198           ENDDO ! of DO iq=1,nq
1199
1200           ! add condensation tendency for H2O2
1201           if (igcm_h2o2.ne.0) then
1202             DO ig=1,ngrid
1203               dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
1204     &                                +zdqscloud(ig,igcm_h2o2)
1205             ENDDO
1206           endif
1207
1208         END IF  ! of IF (photochem.or.thermochem)
1209#endif
1210
1211c   6d. Updates
1212c     ---------
1213
1214        DO iq=1, nq
1215          DO ig=1,ngrid
1216
1217c       ---------------------------------
1218c       Updating tracer budget on surface
1219c       ---------------------------------
1220            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1221
1222          ENDDO  ! (ig)
1223        ENDDO    ! (iq)
1224
1225      endif !  of if (tracer)
1226
1227#ifndef MESOSCALE
1228c-----------------------------------------------------------------------
1229c   7. THERMOSPHERE CALCULATION
1230c-----------------------------------------------------------------------
1231
1232      if (callthermos) then
1233        call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
1234     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1235     &     pt,pq,pu,pv,pdt,pdq,
1236     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1237
1238        DO l=1,nlayer
1239          DO ig=1,ngrid
1240            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
1241            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
1242     &                         +zdteuv(ig,l)
1243            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
1244            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
1245            DO iq=1, nq
1246              pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
1247            ENDDO
1248          ENDDO
1249        ENDDO
1250
1251      endif ! of if (callthermos)
1252#endif
1253
1254c-----------------------------------------------------------------------
1255c   8. Carbon dioxide condensation-sublimation:
1256c     (should be the last atmospherical physical process to be computed)
1257c   -------------------------------------------
1258
1259      IF (tituscap) THEN
1260         !!! get the actual co2 seasonal cap from Titus observations
1261         CALL geticecover( ngrid, 180.*zls/pi,
1262     .                  180.*long/pi, 180.*lati/pi, co2ice )
1263         co2ice = co2ice * 10000.
1264      ENDIF
1265     
1266     
1267      pdpsrf(:) = 0
1268
1269      IF (callcond) THEN
1270         CALL newcondens(ngrid,nlayer,nq,ptimestep,
1271     $              capcal,zplay,zplev,tsurf,pt,
1272     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
1273     $              co2ice,albedo,emis,
1274     $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
1275     $              fluxsurf_sw,zls)
1276
1277         DO l=1,nlayer
1278           DO ig=1,ngrid
1279             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1280             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1281             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1282           ENDDO
1283         ENDDO
1284         DO ig=1,ngrid
1285            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1286         ENDDO
1287
1288         IF (tracer) THEN
1289           DO iq=1, nq
1290            DO l=1,nlayer
1291              DO ig=1,ngrid
1292                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1293              ENDDO
1294            ENDDO
1295           ENDDO
1296         ENDIF ! of IF (tracer)
1297
1298#ifndef MESOSCALE
1299        ! update surface pressure
1300        DO ig=1,ngrid
1301          ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
1302        ENDDO
1303       
1304        ! update pressure levels
1305        DO l=1,nlayer
1306         DO ig=1,ngrid
1307          zplay(ig,l) = aps(l) + bps(l)*ps(ig)
1308          zplev(ig,l) = ap(l)  + bp(l)*ps(ig)
1309         ENDDO
1310        ENDDO
1311        zplev(:,nlayer+1) = 0.
1312       
1313        ! update layers altitude
1314        DO l=2,nlayer
1315          DO ig=1,ngrid
1316           z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
1317           z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
1318           zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
1319          ENDDO
1320        ENDDO
1321#endif
1322     
1323      ENDIF  ! of IF (callcond)
1324     
1325
1326c-----------------------------------------------------------------------
1327c   9. Surface  and sub-surface soil temperature
1328c-----------------------------------------------------------------------
1329c
1330c
1331c   9.1 Increment Surface temperature:
1332c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1333
1334      DO ig=1,ngrid
1335         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1336      ENDDO
1337
1338c  Prescribe a cold trap at south pole (except at high obliquity !!)
1339c  Temperature at the surface is set there to be the temperature
1340c  corresponding to equilibrium temperature between phases of CO2
1341
1342
1343      IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
1344!#ifndef MESOSCALE
1345!         if (caps.and.(obliquit.lt.27.)) then => now done in newcondens
1346           ! NB: Updated surface pressure, at grid point 'ngrid', is
1347           !     ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
1348!           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
1349!     &                     (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
1350!           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*ps(ngrid)))
1351!         endif
1352!#endif
1353c       -------------------------------------------------------------
1354c       Change of surface albedo in case of ground frost
1355c       everywhere except on the north permanent cap and in regions
1356c       covered by dry ice.
1357c              ALWAYS PLACE these lines after newcondens !!!
1358c       -------------------------------------------------------------
1359         do ig=1,ngrid
1360           if ((co2ice(ig).eq.0).and.
1361     &        (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
1362              albedo(ig,1) = albedo_h2o_ice
1363              albedo(ig,2) = albedo_h2o_ice
1364c              write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice)
1365c              write(*,*) "physiq.F frost :"
1366c     &        ,lati(ig)*180./pi, long(ig)*180./pi
1367           endif
1368         enddo  ! of do ig=1,ngrid
1369      ENDIF  ! of IF (tracer.AND.water.AND.(ngrid.NE.1))
1370
1371c
1372c   9.2 Compute soil temperatures and subsurface heat flux:
1373c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1374      IF (callsoil) THEN
1375c       Thermal inertia feedback
1376        IF (tifeedback) THEN
1377         CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
1378         CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
1379     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
1380        ELSE
1381         CALL soil(ngrid,nsoilmx,.false.,inertiedat,
1382     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
1383        ENDIF
1384      ENDIF
1385     
1386
1387c-----------------------------------------------------------------------
1388c  10. Write output files
1389c  ----------------------
1390
1391c    -------------------------------
1392c    Dynamical fields incrementation
1393c    -------------------------------
1394c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
1395      ! temperature, zonal and meridional wind
1396      DO l=1,nlayer
1397        DO ig=1,ngrid
1398          zt(ig,l)=pt(ig,l)  + pdt(ig,l)*ptimestep
1399          zu(ig,l)=pu(ig,l)  + pdu(ig,l)*ptimestep
1400          zv(ig,l)=pv(ig,l)  + pdv(ig,l)*ptimestep
1401        ENDDO
1402      ENDDO
1403
1404      ! tracers
1405      DO iq=1, nq
1406        DO l=1,nlayer
1407          DO ig=1,ngrid
1408            zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
1409          ENDDO
1410        ENDDO
1411      ENDDO
1412
1413      ! Density
1414      DO l=1,nlayer
1415         DO ig=1,ngrid
1416            rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
1417         ENDDO
1418      ENDDO
1419
1420      ! Potential Temperature
1421
1422       DO ig=1,ngrid
1423          DO l=1,nlayer
1424              zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp
1425          ENDDO
1426       ENDDO
1427
1428
1429c    Compute surface stress : (NB: z0 is a common in surfdat.h)
1430c     DO ig=1,ngrid
1431c        cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
1432c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
1433c     ENDDO
1434
1435c     Sum of fluxes in solar spectral bands (for output only)
1436      DO ig=1,ngrid
1437             fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
1438             fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
1439      ENDDO
1440c ******* TEST ******************************************************
1441      ztim1 = 999
1442      DO l=1,nlayer
1443        DO ig=1,ngrid
1444           if (pt(ig,l).lt.ztim1) then
1445               ztim1 = pt(ig,l)
1446               igmin = ig
1447               lmin = l
1448           end if
1449        ENDDO
1450      ENDDO
1451      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
1452        write(*,*) 'PHYSIQ: stability WARNING :'
1453        write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
1454     &              'ig l =', igmin, lmin
1455      end if
1456c *******************************************************************
1457
1458c     ---------------------
1459c     Outputs to the screen
1460c     ---------------------
1461
1462      IF (lwrite) THEN
1463         PRINT*,'Global diagnostics for the physics'
1464         PRINT*,'Variables and their increments x and dx/dt * dt'
1465         WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
1466         WRITE(*,'(2f10.5,2f15.5)')
1467     s   tsurf(igout),zdtsurf(igout)*ptimestep,
1468     s   zplev(igout,1),pdpsrf(igout)*ptimestep
1469         WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
1470         WRITE(*,'(i4,6f10.5)') (l,
1471     s   pu(igout,l),pdu(igout,l)*ptimestep,
1472     s   pv(igout,l),pdv(igout,l)*ptimestep,
1473     s   pt(igout,l),pdt(igout,l)*ptimestep,
1474     s   l=1,nlayer)
1475      ENDIF ! of IF (lwrite)
1476
1477c        ----------------------------------------------------------
1478c        ----------------------------------------------------------
1479c        INTERPOLATIONS IN THE SURFACE-LAYER
1480c        ----------------------------------------------------------
1481c        ----------------------------------------------------------
1482
1483           n_out=0 ! number of elements in the z_out array.
1484                   ! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
1485                   ! to 5
1486           IF (n_out .ne. 0) THEN
1487
1488           IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
1489           IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
1490           IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))
1491
1492           z_out(:)=[3.,2.,1.,0.5,0.1]
1493           u_out(:,:)=0.
1494           T_out(:,:)=0.
1495
1496           call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
1497     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
1498     & T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
1499                   ! pourquoi ustar recalcule ici? fait dans vdifc.
1500
1501#ifndef MESOSCALE
1502            IF (ngrid .eq. 1) THEN
1503               dimout=0
1504            ELSE
1505               dimout=2
1506            ENDIF
1507            DO n=1,n_out
1508               write(zstring, '(F8.6)') z_out(n)
1509               call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
1510     &   'potential temperature at z_out','K',dimout,T_out(:,n))
1511               call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
1512     &   'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
1513            ENDDO
1514            call WRITEDIAGFI(ngrid,'u_star',
1515     &   'friction velocity','m/s',dimout,ustar)
1516            call WRITEDIAGFI(ngrid,'teta_star',
1517     &   'friction potential temperature','K',dimout,tstar)
1518            call WRITEDIAGFI(ngrid,'L',
1519     &   'Monin Obukhov length','m',dimout,L_mo)
1520            call WRITEDIAGFI(ngrid,'vvv',
1521     &   'Vertical velocity variance at zout','m',dimout,vvv)
1522            call WRITEDIAGFI(ngrid,'vhf',
1523     &   'Vertical heat flux at zout','m',dimout,vhf)
1524#else
1525         T_out1(:)=T_out(:,1)
1526         u_out1(:)=u_out(:,1)
1527#endif
1528
1529         ENDIF
1530
1531c        ----------------------------------------------------------
1532c        ----------------------------------------------------------
1533c        END OF SURFACE LAYER INTERPOLATIONS
1534c        ----------------------------------------------------------
1535c        ----------------------------------------------------------
1536
1537      IF (ngrid.NE.1) THEN
1538
1539#ifndef MESOSCALE
1540c        -------------------------------------------------------------------
1541c        Writing NetCDF file  "RESTARTFI" at the end of the run
1542c        -------------------------------------------------------------------
1543c        Note: 'restartfi' is stored just before dynamics are stored
1544c              in 'restart'. Between now and the writting of 'restart',
1545c              there will have been the itau=itau+1 instruction and
1546c              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1547c              thus we store for time=time+dtvr
1548
1549         IF(   ((ecritstart.GT.0) .and.
1550     .          (MOD(icount*iphysiq,ecritstart).EQ.0))
1551     .    .or. lastcall  ) THEN
1552         
1553          ztime_fin = pday + ptime  + ptimestep/(float(iphysiq)*daysec)
1554     .               - day_ini - time_phys
1555          print*, pday,ptime,day_ini, time_phys
1556          write(*,'(A,I7,A,F12.5)')
1557     .         'PHYSIQ: Ecriture du fichier restartfi ; icount=',
1558     .          icount,' date=',ztime_fin
1559           
1560           
1561          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
1562     .                ptimestep,ztime_fin,
1563     .                tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
1564         
1565         ENDIF
1566#endif
1567
1568c        -------------------------------------------------------------------
1569c        Calculation of diagnostic variables written in both stats and
1570c          diagfi files
1571c        -------------------------------------------------------------------
1572
1573         if (tracer) then
1574         
1575           if(doubleq) then
1576              do ig=1,ngrid
1577                dqdustsurf(ig) =
1578     &                zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
1579                dndustsurf(ig) =
1580     &                zdqssed(ig,igcm_dust_number)*tauscaling(ig)
1581                ndust(ig,:) =
1582     &                pq(ig,:,igcm_dust_number)*tauscaling(ig)
1583                qdust(ig,:) =
1584     &                pq(ig,:,igcm_dust_mass)*tauscaling(ig)
1585              enddo
1586              if (scavenging) then
1587                do ig=1,ngrid
1588                  dqdustsurf(ig) = dqdustsurf(ig) +
1589     &                     zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
1590                  dndustsurf(ig) = dndustsurf(ig) +
1591     &                     zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
1592                  nccn(ig,:) =
1593     &                     pq(ig,:,igcm_ccn_number)*tauscaling(ig)
1594                  qccn(ig,:) =
1595     &                     pq(ig,:,igcm_ccn_mass)*tauscaling(ig)
1596                enddo
1597              endif
1598           endif
1599                   
1600           if (water) then
1601             mtot(:)=0
1602             icetot(:)=0
1603             rave(:)=0
1604             tauTES(:)=0
1605             do ig=1,ngrid
1606               do l=1,nlayer
1607                 mtot(ig) = mtot(ig) +
1608     &                      zq(ig,l,igcm_h2o_vap) *
1609     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
1610                 icetot(ig) = icetot(ig) +
1611     &                        zq(ig,l,igcm_h2o_ice) *
1612     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
1613c                Computing abs optical depth at 825 cm-1 in each
1614c                  layer to simulate NEW TES retrieval
1615                 Qabsice = min(
1616     &             max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2
1617     &                        )
1618                 opTES(ig,l)= 0.75 * Qabsice *
1619     &             zq(ig,l,igcm_h2o_ice) *
1620     &             (zplev(ig,l) - zplev(ig,l+1)) / g
1621     &             / (rho_ice * rice(ig,l) * (1.+nuice_ref))
1622                 tauTES(ig)=tauTES(ig)+ opTES(ig,l)
1623               enddo
1624c              rave(ig)=rave(ig)/max(icetot(ig),1.e-30)       ! mass weight
1625c               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1626             enddo
1627             call watersat(ngrid*nlayer,zt,zplay,zqsat)
1628             satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
1629
1630             if (scavenging) then
1631               Nccntot(:)= 0
1632               Mccntot(:)= 0
1633               rave(:)=0
1634               do ig=1,ngrid
1635                 do l=1,nlayer
1636                    Nccntot(ig) = Nccntot(ig) +
1637     &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
1638     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
1639                    Mccntot(ig) = Mccntot(ig) +
1640     &              zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
1641     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
1642cccc Column integrated effective ice radius
1643cccc is weighted by total ice surface area (BETTER than total ice mass)
1644                    rave(ig) = rave(ig) +
1645     &                      tauscaling(ig) *
1646     &                      zq(ig,l,igcm_ccn_number) *
1647     &                      (zplev(ig,l) - zplev(ig,l+1)) / g *
1648     &                      rice(ig,l) * rice(ig,l)*  (1.+nuice_ref)
1649                 enddo
1650               rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
1651     &                  /max(pi*rave(ig),1.e-30) ! surface weight
1652               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
1653               enddo
1654             else ! of if (scavenging)
1655               rave(:)=0
1656               do ig=1,ngrid
1657                 do l=1,nlayer
1658                 rave(ig) = rave(ig) +
1659     &                      zq(ig,l,igcm_h2o_ice) *
1660     &                      (zplev(ig,l) - zplev(ig,l+1)) / g *
1661     &                      rice(ig,l) * (1.+nuice_ref)
1662                 enddo
1663                 rave(ig) = max(rave(ig) /
1664     &             max(icetot(ig),1.e-30),1.e-30) ! mass weight
1665               enddo
1666             endif ! of if (scavenging)
1667
1668           endif ! of if (water)
1669         endif ! of if (tracer)
1670
1671#ifndef MESOSCALE
1672c        -----------------------------------------------------------------
1673c        WSTATS: Saving statistics
1674c        -----------------------------------------------------------------
1675c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
1676c        which can later be used to make the statistic files of the run:
1677c        "stats")          only possible in 3D runs !
1678         
1679       IF (callstats) THEN
1680
1681        call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1682        call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1683        call wstats(ngrid,"co2ice","CO2 ice cover",
1684     &                "kg.m-2",2,co2ice)
1685        call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
1686     &                2,tauref)
1687        call wstats(ngrid,"fluxsurf_lw",
1688     &                "Thermal IR radiative flux to surface","W.m-2",2,
1689     &                fluxsurf_lw)
1690        call wstats(ngrid,"fluxsurf_sw",
1691     &                "Solar radiative flux to surface","W.m-2",2,
1692     &                fluxsurf_sw_tot)
1693        call wstats(ngrid,"fluxtop_lw",
1694     &                "Thermal IR radiative flux to space","W.m-2",2,
1695     &                fluxtop_lw)
1696        call wstats(ngrid,"fluxtop_sw",
1697     &                "Solar radiative flux to space","W.m-2",2,
1698     &                fluxtop_sw_tot)
1699        call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1700        call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1701        call wstats(ngrid,"v","Meridional (North-South) wind",
1702     &                "m.s-1",3,zv)
1703        call wstats(ngrid,"w","Vertical (down-up) wind",
1704     &                "m.s-1",3,pw)
1705        call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
1706        call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
1707          call wstats(ngrid,"q2",
1708     &                "Boundary layer eddy kinetic energy",
1709     &                "m2.s-2",3,q2)
1710          call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
1711     &                emis)
1712c          call wstats(ngrid,"ssurf","Surface stress","N.m-2",
1713c    &                2,zstress)
1714c          call wstats(ngrid,"sw_htrt","sw heat.rate",
1715c    &                 "W.m-2",3,zdtsw)
1716c          call wstats(ngrid,"lw_htrt","lw heat.rate",
1717c    &                 "W.m-2",3,zdtlw)
1718
1719          if (calltherm) then
1720            call wstats(ngrid,"zmax_th","Height of thermals",
1721     &                "m",2,zmax_th)
1722            call wstats(ngrid,"hfmax_th","Max thermals heat flux",
1723     &                "K.m/s",2,hfmax_th)
1724            call wstats(ngrid,"wstar",
1725     &                "Max vertical velocity in thermals",
1726     &                "m/s",2,wstar)
1727          endif
1728
1729           if (tracer) then
1730             if (water) then
1731               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
1732     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
1733               call wstats(ngrid,"vmr_h2ovap",
1734     &                    "H2O vapor volume mixing ratio","mol/mol",
1735     &                    3,vmr)
1736               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1737     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
1738               call wstats(ngrid,"vmr_h2oice",
1739     &                    "H2O ice volume mixing ratio","mol/mol",
1740     &                    3,vmr)
1741              ! also store vmr_ice*rice for better diagnostics of rice
1742               vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
1743     &                               zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1744               call wstats(ngrid,"vmr_h2oice_rice",
1745     &                "H2O ice mixing ratio times ice particule size",
1746     &                    "(mol/mol)*m",
1747     &                    3,vmr)
1748               vmr=zqsat(1:ngrid,1:nlayer)
1749     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
1750               call wstats(ngrid,"vmr_h2osat",
1751     &                    "saturation volume mixing ratio","mol/mol",
1752     &                    3,vmr)
1753               call wstats(ngrid,"h2o_ice_s",
1754     &                    "surface h2o_ice","kg/m2",
1755     &                    2,qsurf(1,igcm_h2o_ice))
1756               call wstats(ngrid,'albedo',
1757     &                         'albedo',
1758     &                         '',2,albedo(1,1))
1759               call wstats(ngrid,"mtot",
1760     &                    "total mass of water vapor","kg/m2",
1761     &                    2,mtot)
1762               call wstats(ngrid,"icetot",
1763     &                    "total mass of water ice","kg/m2",
1764     &                    2,icetot)
1765               call wstats(ngrid,"reffice",
1766     &                    "Mean reff","m",
1767     &                    2,rave)
1768              call wstats(ngrid,"Nccntot",
1769     &                    "condensation nuclei","Nbr/m2",
1770     &                    2,Nccntot)
1771              call wstats(ngrid,"Mccntot",
1772     &                    "condensation nuclei mass","kg/m2",
1773     &                    2,Mccntot)
1774              call wstats(ngrid,"rice",
1775     &                    "Ice particle size","m",
1776     &                    3,rice)
1777               if (.not.activice) then
1778                 call wstats(ngrid,"tauTESap",
1779     &                      "tau abs 825 cm-1","",
1780     &                      2,tauTES)
1781               else
1782                 call wstats(ngrid,'tauTES',
1783     &                   'tau abs 825 cm-1',
1784     &                  '',2,taucloudtes)
1785               endif
1786
1787             endif ! of if (water)
1788             
1789             
1790           if (dustbin.ne.0) then
1791         
1792             call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
1793             
1794             if (doubleq) then
1795c            call wstats(ngrid,'qsurf','qsurf',
1796c     &                       'kg.m-2',2,qsurf(1,igcm_dust_mass))
1797c            call wstats(ngrid,'Nsurf','N particles',
1798c     &                       'N.m-2',2,qsurf(1,igcm_dust_number))
1799c            call wstats(ngrid,'dqsdev','ddevil lift',
1800c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
1801c            call wstats(ngrid,'dqssed','sedimentation',
1802c     &                       'kg.m-2.s-1',2,zdqssed(1,1))
1803c            call wstats(ngrid,'dqsdif','diffusion',
1804c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
1805               call wstats(ngrid,'dqsdust',
1806     &                        'deposited surface dust mass',
1807     &                        'kg.m-2.s-1',2,dqdustsurf)
1808               call wstats(ngrid,'dqndust',
1809     &                        'deposited surface dust number',
1810     &                        'number.m-2.s-1',2,dndustsurf)
1811               call wstats(ngrid,'reffdust','reffdust',
1812     &                        'm',3,rdust*ref_r0)
1813               call wstats(ngrid,'dustq','Dust mass mr',
1814     &                        'kg/kg',3,qdust)
1815               call wstats(ngrid,'dustN','Dust number',
1816     &                        'part/kg',3,ndust)
1817             else
1818               do iq=1,dustbin
1819                 write(str2(1:2),'(i2.2)') iq
1820                 call wstats(ngrid,'q'//str2,'mix. ratio',
1821     &                         'kg/kg',3,zq(1,1,iq))
1822                 call wstats(ngrid,'qsurf'//str2,'qsurf',
1823     &                         'kg.m-2',2,qsurf(1,iq))
1824               end do
1825             endif ! (doubleq)
1826
1827             if (scavenging) then
1828               call wstats(ngrid,'ccnq','CCN mass mr',
1829     &                        'kg/kg',3,qccn)
1830               call wstats(ngrid,'ccnN','CCN number',
1831     &                        'part/kg',3,nccn)
1832             endif ! (scavenging)
1833         
1834          endif ! (dustbin.ne.0)
1835
1836
1837             if (thermochem.or.photochem) then
1838                do iq=1,nq
1839                  if (noms(iq) .ne. "dust_mass" .and.
1840     $                 noms(iq) .ne. "dust_number" .and.
1841     $                 noms(iq) .ne. "ccn_mass" .and.
1842     $                 noms(iq) .ne. "ccn_number" .and.
1843     $                 noms(iq) .ne. "h2o_vap" .and.
1844     $                 noms(iq) .ne. "h2o_ice") then
1845                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
1846     &                          *mmean(1:ngrid,1:nlayer)/mmol(iq)
1847                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
1848     &                          *rho(1:ngrid,1:nlayer)*n_avog/
1849     &                           (1000*mmol(iq))
1850                   call wstats(ngrid,"vmr_"//trim(noms(iq)),
1851     $                         "Volume mixing ratio","mol/mol",3,vmr)
1852!                   call wstats(ngrid,"rho_"//trim(noms(iq)),
1853!     $                     "Number density","cm-3",3,rhopart)
1854!                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
1855!     $                     "Number density","cm-3",3,rhopart)
1856                   do ig = 1,ngrid
1857                      colden(ig,iq) = 0.                           
1858                   end do
1859                   do l=1,nlayer                                   
1860                      do ig=1,ngrid                                 
1861                         colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
1862     $                                  *(zplev(ig,l)-zplev(ig,l+1))
1863     $                                  *6.022e22/(mmol(iq)*g)       
1864                      end do                                       
1865                   end do                                         
1866                   call wstats(ngrid,"c_"//trim(noms(iq)),           
1867     $                         "column","mol cm-2",2,colden(1,iq)) 
1868                  end if ! of if (noms(iq) .ne. "dust_mass" ...)
1869                end do ! of do iq=1,nq
1870             end if ! of if (thermochem.or.photochem)
1871
1872           end if ! of if (tracer)
1873
1874           IF(lastcall) THEN
1875             write (*,*) "Writing stats..."
1876             call mkstats(ierr)
1877           ENDIF
1878
1879         ENDIF !if callstats
1880
1881c        (Store EOF for Mars Climate database software)
1882         IF (calleofdump) THEN
1883            CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
1884         ENDIF
1885#endif
1886!endif of ifndef MESOSCALE
1887
1888#ifdef MESOSCALE
1889     
1890      !! see comm_wrf.
1891      !! not needed when an array is already in a shared module.
1892      !! --> example : hfmax_th, zmax_th
1893
1894      !state  real  HR_SW     ikj   misc  1  -  h  "HR_SW"     "HEATING RATE SW"                 "K/s"
1895      comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
1896      !state  real  HR_LW     ikj   misc  1  -  h  "HR_LW"     "HEATING RATE LW"                 "K/s"
1897      comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
1898      !state  real  SWDOWNZ    ij   misc  1  -  h  "SWDOWNZ"   "DOWNWARD SW FLUX AT SURFACE"     "W m-2"
1899      comm_SWDOWNZ(1:ngrid) = fluxsurf_sw_tot(1:ngrid)
1900      !state  real  TAU_DUST   ij   misc  1  -  h  "TAU_DUST"  "REFERENCE VISIBLE DUST OPACITY"  ""
1901      comm_TAU_DUST(1:ngrid) = tauref(1:ngrid)
1902      !state  real  RDUST     ikj   misc  1  -  h  "RDUST"     "DUST RADIUS"                     "m"
1903      comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
1904      !state  real  QSURFDUST  ij   misc  1  -  h  "QSURFDUST" "DUST MASS AT SURFACE"            "kg m-2"
1905      IF (igcm_dust_mass .ne. 0) THEN
1906        comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
1907      ELSE
1908        comm_QSURFDUST(1:ngrid) = 0.
1909      ENDIF
1910      !state  real  MTOT       ij   misc  1  -  h  "MTOT"      "TOTAL MASS WATER VAPOR in pmic"  "pmic"
1911      comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
1912      !state  real  ICETOT     ij   misc  1  -  h  "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"
1913      comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
1914      !state  real  VMR_ICE   ikj   misc  1  -  h  "VMR_ICE"   "VOL. MIXING RATIO ICE"           "ppm"
1915      IF (igcm_h2o_ice .ne. 0) THEN
1916        comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
1917     .    * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
1918     .    * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
1919      ELSE
1920        comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
1921      ENDIF
1922      !state  real  TAU_ICE    ij   misc  1  -  h  "TAU_ICE"   "CLOUD OD at 825 cm-1 TES"        ""
1923      if (activice) then
1924        comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
1925      else
1926        comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
1927      endif
1928      !state  real  RICE      ikj   misc  1  -  h  "RICE"      "ICE RADIUS"                      "m"
1929      comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
1930
1931   
1932      !! calculate sensible heat flux in W/m2 for outputs
1933      !! -- the one computed in vdifc is not the real one
1934      !! -- vdifc must have been called
1935      if (.not.callrichsl) then
1936        sensibFlux(1:ngrid) = zflubid(1:ngrid)
1937     .         - capcal(1:ngrid)*zdtsdif(1:ngrid)
1938      else
1939        sensibFlux(1:ngrid) =
1940     &   (pplay(1:ngrid,1)/(r*pt(1:ngrid,1)))*cpp
1941     &   *sqrt(pu(1:ngrid,1)*pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
1942     &         +(log(1.+0.7*wstar(1:ngrid) + 2.3*wstar(1:ngrid)**2))**2)
1943     &   *zcdh(1:ngrid)*(tsurf(1:ngrid)-zh(1:ngrid,1))
1944      endif
1945         
1946#else
1947#ifndef MESOINI
1948
1949c        ==========================================================
1950c        WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
1951c          any variable for diagnostic (output with period
1952c          "ecritphy", set in "run.def")
1953c        ==========================================================
1954c        WRITEDIAGFI can ALSO be called from any other subroutines
1955c        for any variables !!
1956c        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
1957c    &                  emis)
1958c        call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
1959c        call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
1960         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
1961     &                  tsurf)
1962         call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
1963         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
1964     &                                         ,"kg.m-2",2,co2ice)
1965
1966         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
1967     &                  "K",2,zt(1,7))
1968         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
1969     &                  fluxsurf_lw)
1970         call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
1971     &                  fluxsurf_sw_tot)
1972         call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
1973     &                  fluxtop_lw)
1974         call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
1975     &                  fluxtop_sw_tot)
1976         call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
1977         call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
1978         call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
1979         call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
1980         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
1981c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
1982c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
1983         call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
1984c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
1985c    &                  zstress)
1986c        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
1987c    &                   'w.m-2',3,zdtsw)
1988c        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
1989c    &                   'w.m-2',3,zdtlw)
1990            if (.not.activice) then
1991               CALL WRITEDIAGFI(ngrid,'tauTESap',
1992     &                         'tau abs 825 cm-1',
1993     &                         '',2,tauTES)
1994             else
1995               CALL WRITEDIAGFI(ngrid,'tauTES',
1996     &                         'tau abs 825 cm-1',
1997     &                         '',2,taucloudtes)
1998             endif
1999
2000#else
2001     !!! this is to ensure correct initialisation of mesoscale model
2002        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
2003     &                  tsurf)
2004        call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
2005        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
2006     &                  co2ice)
2007        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
2008        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
2009        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
2010        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
2011     &                  emis)
2012        call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
2013     &                       "K",3,tsoil)
2014        call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
2015     &                       "K",3,inertiedat)
2016#endif
2017
2018
2019c        ----------------------------------------------------------
2020c        Outputs of the CO2 cycle
2021c        ----------------------------------------------------------
2022
2023         if (tracer.and.(igcm_co2.ne.0)) then
2024!          call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
2025!    &                     "kg/kg",2,zq(1,1,igcm_co2))
2026          call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
2027     &                     "kg/kg",3,zq(1,1,igcm_co2))
2028       
2029         ! Compute co2 column
2030         co2col(:)=0
2031         do l=1,nlayer
2032           do ig=1,ngrid
2033             co2col(ig)=co2col(ig)+
2034     &                  zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
2035           enddo
2036         enddo
2037         call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
2038     &                  co2col)
2039         endif ! of if (tracer.and.(igcm_co2.ne.0))
2040
2041c        ----------------------------------------------------------
2042c        Outputs of the water cycle
2043c        ----------------------------------------------------------
2044         if (tracer) then
2045           if (water) then
2046
2047#ifdef MESOINI
2048            !!!! waterice = q01, voir readmeteo.F90
2049            call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
2050     &                      'kg/kg',3,
2051     &                       zq(1:ngrid,1:nlayer,igcm_h2o_ice))
2052            !!!! watervapor = q02, voir readmeteo.F90
2053            call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
2054     &                      'kg/kg',3,
2055     &                       zq(1:ngrid,1:nlayer,igcm_h2o_vap))
2056            !!!! surface waterice qsurf02 (voir readmeteo)
2057            call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
2058     &                      'kg.m-2',2,
2059     &                       qsurf(1:ngrid,igcm_h2o_ice))
2060#endif
2061
2062            CALL WRITEDIAGFI(ngrid,'mtot',
2063     &                       'total mass of water vapor',
2064     &                       'kg/m2',2,mtot)
2065            CALL WRITEDIAGFI(ngrid,'icetot',
2066     &                       'total mass of water ice',
2067     &                       'kg/m2',2,icetot)
2068            vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
2069     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
2070            call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
2071     &                       'mol/mol',3,vmr)
2072            vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
2073     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
2074            call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
2075     &                       'mol/mol',3,vmr)
2076            CALL WRITEDIAGFI(ngrid,'reffice',
2077     &                       'Mean reff',
2078     &                       'm',2,rave)
2079            if (scavenging) then
2080              CALL WRITEDIAGFI(ngrid,"Nccntot",
2081     &                    "condensation nuclei","Nbr/m2",
2082     &                    2,Nccntot)
2083              CALL WRITEDIAGFI(ngrid,"Mccntot",
2084     &                    "mass condensation nuclei","kg/m2",
2085     &                    2,Mccntot)
2086            endif
2087            call WRITEDIAGFI(ngrid,'rice','Ice particle size',
2088     &                       'm',3,rice)
2089            call WRITEDIAGFI(ngrid,'h2o_ice_s',
2090     &                       'surface h2o_ice',
2091     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
2092            CALL WRITEDIAGFI(ngrid,'albedo',
2093     &                         'albedo',
2094     &                         '',2,albedo(1,1))
2095              if (tifeedback) then
2096                 call WRITEDIAGSOIL(ngrid,"soiltemp",
2097     &                              "Soil temperature","K",
2098     &                              3,tsoil)
2099                 call WRITEDIAGSOIL(ngrid,'soilti',
2100     &                       'Soil Thermal Inertia',
2101     &                       'J.s-1/2.m-2.K-1',3,inertiesoil)
2102              endif
2103           endif !(water)
2104
2105
2106           if (water.and..not.photochem) then
2107             iq=nq
2108c            write(str2(1:2),'(i2.2)') iq
2109c            call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
2110c    &                       'kg.m-2',2,zdqscloud(1,iq))
2111c            call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
2112c    &                       'kg/kg',3,zdqchim(1,1,iq))
2113c            call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
2114c    &                       'kg/kg',3,zdqdif(1,1,iq))
2115c            call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
2116c    &                       'kg/kg',3,zdqadj(1,1,iq))
2117c            call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
2118c    &                       'kg/kg',3,zdqc(1,1,iq))
2119           endif  !(water.and..not.photochem)
2120         endif
2121
2122c        ----------------------------------------------------------
2123c        Outputs of the dust cycle
2124c        ----------------------------------------------------------
2125
2126        call WRITEDIAGFI(ngrid,'tauref',
2127     &                    'Dust ref opt depth','NU',2,tauref)
2128
2129         if (tracer.and.(dustbin.ne.0)) then
2130
2131          call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
2132
2133#ifndef MESOINI
2134           if (doubleq) then
2135c            call WRITEDIAGFI(ngrid,'qsurf','qsurf',
2136c     &                       'kg.m-2',2,qsurf(1,igcm_dust_mass))
2137c            call WRITEDIAGFI(ngrid,'Nsurf','N particles',
2138c     &                       'N.m-2',2,qsurf(1,igcm_dust_number))
2139c            call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
2140c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
2141c            call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
2142c     &                       'kg.m-2.s-1',2,zdqssed(1,1))
2143c            call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
2144c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
2145c             call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
2146c     &                       'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
2147c             call WRITEDIAGFI(ngrid,'subice','sublimated ice',
2148c     &                       'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
2149             call WRITEDIAGFI(ngrid,'dqsdust',
2150     &                        'deposited surface dust mass',
2151     &                        'kg.m-2.s-1',2,dqdustsurf)
2152             call WRITEDIAGFI(ngrid,'dqndust',
2153     &                        'deposited surface dust number',
2154     &                        'number.m-2.s-1',2,dndustsurf)
2155             call WRITEDIAGFI(ngrid,'reffdust','reffdust',
2156     &                        'm',3,rdust*ref_r0)
2157             call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2158     &                        'kg/kg',3,qdust)
2159             call WRITEDIAGFI(ngrid,'dustN','Dust number',
2160     &                        'part/kg',3,ndust)
2161c             call WRITEDIAGFI(ngrid,"tauscaling",
2162c     &                    "dust conversion factor"," ",2,tauscaling)
2163           else
2164             do iq=1,dustbin
2165               write(str2(1:2),'(i2.2)') iq
2166               call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio',
2167     &                         'kg/kg',3,zq(1,1,iq))
2168               call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
2169     &                         'kg.m-2',2,qsurf(1,iq))
2170             end do
2171           endif ! (doubleq)
2172
2173           if (scavenging) then
2174             call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
2175     &                        'kg/kg',3,qccn)
2176             call WRITEDIAGFI(ngrid,'ccnN','CCN number',
2177     &                        'part/kg',3,nccn)
2178           endif ! (scavenging)
2179
2180c          if (submicron) then
2181c            call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
2182c    &                        'kg/kg',3,pq(1,1,igcm_dust_submicron))
2183c          endif ! (submicron)
2184
2185#else
2186!     !!! to initialize mesoscale we need scaled variables
2187!     !!! because this must correspond to starting point for tracers
2188!             call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2189!     &           'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
2190!             call WRITEDIAGFI(ngrid,'dustN','Dust number',
2191!     &           'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
2192!             call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
2193!     &           'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
2194!             call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
2195!     &           'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
2196      if (freedust) then
2197             call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2198     &                        'kg/kg',3,qdust)
2199             call WRITEDIAGFI(ngrid,'dustN','Dust number',
2200     &                        'part/kg',3,ndust)
2201             call WRITEDIAGFI(ngrid,'ccn','CCN mass mr',
2202     &                        'kg/kg',3,qccn)
2203             call WRITEDIAGFI(ngrid,'ccnN','CCN number',
2204     &                        'part/kg',3,nccn)
2205      else
2206             call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
2207     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
2208             call WRITEDIAGFI(ngrid,'dustN','Dust number',
2209     &                        'part/kg',3,pq(1,1,igcm_dust_number))
2210             call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
2211     &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
2212             call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
2213     &                        'part/kg',3,pq(1,1,igcm_ccn_number))
2214      endif
2215#endif
2216
2217         end if  ! (tracer.and.(dustbin.ne.0))
2218
2219
2220c        ----------------------------------------------------------
2221c        Thermospheric outputs
2222c        ----------------------------------------------------------
2223
2224         if(callthermos) then
2225
2226            call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
2227     $           3,zdtnlte)
2228            call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
2229     $           3,zdteuv)
2230            call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
2231     $           3,zdtconduc)
2232            call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
2233     $           3,zdtnirco2)
2234
2235
2236             if (thermochem.or.photochem) then
2237                do iq=1,nq
2238                  if (noms(iq) .ne. "dust_mass" .and.
2239     $                 noms(iq) .ne. "dust_number" .and.
2240     $                 noms(iq) .ne. "ccn_mass" .and.
2241     $                 noms(iq) .ne. "ccn_number" .and.
2242     $                 noms(iq) .ne. "h2o_vap" .and.
2243     $                 noms(iq) .ne. "h2o_ice") then
2244                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
2245     &                          *mmean(1:ngrid,1:nlayer)/mmol(iq)
2246                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
2247     &                          *rho(1:ngrid,1:nlayer)*n_avog/
2248     &                           (1000*mmol(iq))
2249                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
2250     $                     "Number density","cm-3",3,rhopart)
2251                   if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
2252     $                 (noms(iq).eq."o3")) then
2253                      call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
2254     $                         "Volume mixing ratio","mol/mol",3,vmr)
2255                   end if
2256                   do ig = 1,ngrid
2257                      colden(ig,iq) = 0.
2258                   end do
2259                   do l=1,nlayer
2260                      do ig=1,ngrid
2261                         colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
2262     $                                  *(zplev(ig,l)-zplev(ig,l+1))
2263     $                                  *6.022e22/(mmol(iq)*g)
2264                      end do
2265                   end do
2266                   call writediagfi(ngrid,"c_"//trim(noms(iq)),
2267     $                             "column","mol cm-2",2,colden(1,iq))
2268                  end if ! of if (noms(iq) .ne. "dust_mass" ...)
2269                end do ! of do iq=1,nq
2270             end if ! of if (thermochem.or.photochem)
2271
2272
2273         endif  !(callthermos)
2274
2275c        ----------------------------------------------------------
2276c        ----------------------------------------------------------
2277c        PBL OUTPUS
2278c        ----------------------------------------------------------
2279c        ----------------------------------------------------------
2280
2281c        ----------------------------------------------------------
2282c        Outputs of thermals
2283c        ----------------------------------------------------------
2284         if (calltherm) then
2285
2286!        call WRITEDIAGFI(ngrid,'dtke',
2287!     &              'tendance tke thermiques','m**2/s**2',
2288!     &                         3,dtke_th)
2289!        call WRITEDIAGFI(ngrid,'d_u_ajs',
2290!     &              'tendance u thermiques','m/s',
2291!     &                         3,pdu_th*ptimestep)
2292!        call WRITEDIAGFI(ngrid,'d_v_ajs',
2293!     &              'tendance v thermiques','m/s',
2294!     &                         3,pdv_th*ptimestep)
2295!        if (tracer) then
2296!        if (nq .eq. 2) then
2297!        call WRITEDIAGFI(ngrid,'deltaq_th',
2298!     &              'delta q thermiques','kg/kg',
2299!     &                         3,ptimestep*pdq_th(:,:,2))
2300!        endif
2301!        endif
2302
2303        call WRITEDIAGFI(ngrid,'zmax_th',
2304     &              'hauteur du thermique','m',
2305     &                         2,zmax_th)
2306        call WRITEDIAGFI(ngrid,'hfmax_th',
2307     &              'maximum TH heat flux','K.m/s',
2308     &                         2,hfmax_th)
2309        call WRITEDIAGFI(ngrid,'wstar',
2310     &              'maximum TH vertical velocity','m/s',
2311     &                         2,wstar)
2312
2313         endif
2314
2315c        ----------------------------------------------------------
2316c        ----------------------------------------------------------
2317c        END OF PBL OUTPUS
2318c        ----------------------------------------------------------
2319c        ----------------------------------------------------------
2320
2321
2322c        ----------------------------------------------------------
2323c        Output in netcdf file "diagsoil.nc" for subterranean
2324c          variables (output every "ecritphy", as for writediagfi)
2325c        ----------------------------------------------------------
2326
2327         ! Write soil temperature
2328!        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
2329!    &                     3,tsoil)
2330         ! Write surface temperature
2331!        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
2332!    &                     2,tsurf)
2333
2334c        ==========================================================
2335c        END OF WRITEDIAGFI
2336c        ==========================================================
2337#endif
2338! of ifdef MESOSCALE
2339
2340      ELSE     ! if(ngrid.eq.1)
2341
2342#ifndef MESOSCALE
2343         write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
2344     &    zls*180./pi,odpref,tauref
2345c      ----------------------------------------------------------------------
2346c      Output in grads file "g1d" (ONLY when using testphys1d)
2347c      (output at every X physical timestep)
2348c      ----------------------------------------------------------------------
2349c
2350c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
2351c         CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
2352c         CALL writeg1d(ngrid,1,ps,'ps','Pa')
2353         
2354c         CALL writeg1d(ngrid,nlayer,zt,'T','K')
2355c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
2356c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
2357c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
2358
2359! THERMALS STUFF 1D
2360         if(calltherm) then
2361
2362        call WRITEDIAGFI(ngrid,'lmax_th',
2363     &              'hauteur du thermique','point',
2364     &                         0,lmax_th_out)
2365        call WRITEDIAGFI(ngrid,'zmax_th',
2366     &              'hauteur du thermique','m',
2367     &                         0,zmax_th)
2368        call WRITEDIAGFI(ngrid,'hfmax_th',
2369     &              'maximum TH heat flux','K.m/s',
2370     &                         0,hfmax_th)
2371        call WRITEDIAGFI(ngrid,'wstar',
2372     &              'maximum TH vertical velocity','m/s',
2373     &                         0,wstar)
2374
2375         co2col(:)=0.
2376         if (tracer) then
2377         do l=1,nlayer
2378           do ig=1,ngrid
2379             co2col(ig)=co2col(ig)+
2380     &                  zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
2381         enddo
2382         enddo
2383
2384         end if
2385         call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass'          &
2386     &                                      ,'kg/m-2',0,co2col)
2387         endif ! of if (calltherm)
2388
2389         call WRITEDIAGFI(ngrid,'w','vertical velocity'                 &
2390     &                              ,'m/s',1,pw)
2391         call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
2392         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
2393     &                  tsurf)
2394         call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
2395         call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
2396
2397         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
2398         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
2399         call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
2400!         call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",              &
2401!     &              "K.s-1",1,dtrad/zpopsk)
2402!        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
2403!     &                   'w.m-2',1,zdtsw/zpopsk)
2404!        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
2405!     &                   'w.m-2',1,zdtlw/zpopsk)
2406         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
2407     &                                   ,"kg.m-2",0,co2ice)
2408
2409! or output in diagfi.nc (for testphys1d)
2410         call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
2411         call WRITEDIAGFI(ngrid,'temp','Temperature',
2412     &                       'K',1,zt)
2413     
2414         if(tracer) then
2415c           CALL writeg1d(ngrid,1,tau,'tau','SI')
2416            do iq=1,nq
2417c              CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
2418               call WRITEDIAGFI(ngrid,trim(noms(iq)),
2419     &              trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
2420            end do
2421           if (doubleq) then
2422             call WRITEDIAGFI(ngrid,'rdust','rdust',
2423     &                        'm',1,rdust)
2424           endif
2425           if (water.AND.tifeedback) then
2426             call WRITEDIAGFI(ngrid,"soiltemp",
2427     &                              "Soil temperature","K",
2428     &                              1,tsoil)
2429             call WRITEDIAGFI(ngrid,'soilti',
2430     &                       'Soil Thermal Inertia',
2431     &                       'J.s-1/2.m-2.K-1',1,inertiesoil)
2432           endif
2433         end if
2434         
2435cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
2436ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2437         IF (water) THEN
2438
2439           if (.not.activice) then
2440
2441             tauTES=0
2442             do l=1,nlayer
2443               Qabsice = min(
2444     &             max(0.4e6*rice(1,l)*(1.+nuice_ref)-0.05 ,0.),1.2
2445     &                      )
2446               opTES(1,l)= 0.75 * Qabsice *
2447     &             zq(1,l,igcm_h2o_ice) *
2448     &             (zplev(1,l) - zplev(1,l+1)) / g
2449     &             / (rho_ice * rice(1,l) * (1.+nuice_ref))
2450               tauTES=tauTES+ opTES(1,l)
2451             enddo
2452             CALL WRITEDIAGFI(ngrid,'tauTESap',
2453     &                         'tau abs 825 cm-1',
2454     &                         '',0,tauTES)
2455           else
2456
2457             CALL WRITEDIAGFI(ngrid,'tauTES',
2458     &                         'tau abs 825 cm-1',
2459     &                         '',0,taucloudtes)
2460           endif
2461     
2462           mtot = 0
2463           icetot = 0
2464           h2otot = qsurf(1,igcm_h2o_ice)
2465
2466           do l=1,nlayer
2467             mtot = mtot +  zq(1,l,igcm_h2o_vap)
2468     &                 * (zplev(1,l) - zplev(1,l+1)) / g
2469             icetot = icetot +  zq(1,l,igcm_h2o_ice)
2470     &                 * (zplev(1,l) - zplev(1,l+1)) / g
2471           end do
2472           h2otot = h2otot+mtot+icetot
2473
2474             CALL WRITEDIAGFI(ngrid,'h2otot',
2475     &                         'h2otot',
2476     &                         'kg/m2',0,h2otot)
2477             CALL WRITEDIAGFI(ngrid,'mtot',
2478     &                         'mtot',
2479     &                         'kg/m2',0,mtot)
2480             CALL WRITEDIAGFI(ngrid,'icetot',
2481     &                         'icetot',
2482     &                         'kg/m2',0,icetot)
2483
2484           if (scavenging) then
2485
2486             rave = 0
2487             do l=1,nlayer
2488cccc Column integrated effective ice radius
2489cccc is weighted by total ice surface area (BETTER)
2490             rave = rave + tauscaling(1) *
2491     &              zq(1,l,igcm_ccn_number) *
2492     &              (zplev(1,l) - zplev(1,l+1)) / g *
2493     &              rice(1,l) * rice(1,l)*  (1.+nuice_ref)
2494             enddo
2495             rave=icetot*0.75/max(rave*pi*rho_ice,1.e-30) ! surface weight
2496
2497              Nccntot= 0
2498              call watersat(ngrid*nlayer,zt,zplay,zqsat)
2499               do l=1,nlayer
2500                   Nccntot = Nccntot +
2501     &              zq(1,l,igcm_ccn_number)*tauscaling(1)
2502     &              *(zplev(1,l) - zplev(1,l+1)) / g
2503                   satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
2504                   satu(1,l) = (max(satu(1,l),float(1))-1)
2505!     &                      * zq(1,l,igcm_h2o_vap) *
2506!     &                      (zplev(1,l) - zplev(1,l+1)) / g
2507               enddo
2508               call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
2509     &                    satu)
2510               CALL WRITEDIAGFI(ngrid,'Nccntot',
2511     &                         'Nccntot',
2512     &                         'nbr/m2',0,Nccntot)
2513
2514             call WRITEDIAGFI(ngrid,'zdqsed_dustq'
2515     & ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
2516             call WRITEDIAGFI(ngrid,'zdqsed_dustN'
2517     &,'sedimentation N','Nbr.m-2.s-1',1,
2518     &                                 zdqsed(1,:,igcm_dust_number))
2519
2520           else ! of if (scavenging)
2521
2522cccc Column integrated effective ice radius
2523cccc is weighted by total ice mass         (LESS GOOD)
2524             rave = 0
2525             do l=1,nlayer
2526               rave = rave + zq(1,l,igcm_h2o_ice)
2527     &              * (zplev(1,l) - zplev(1,l+1)) / g
2528     &              *  rice(1,l) * (1.+nuice_ref)
2529             enddo
2530             rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
2531           endif ! of if (scavenging)
2532
2533           CALL WRITEDIAGFI(ngrid,'reffice',
2534     &                      'reffice',
2535     &                      'm',0,rave)
2536
2537           do iq=1,nq
2538             call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
2539     &            trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
2540           end do
2541     
2542        call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
2543     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
2544        call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap',
2545     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
2546        call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice',
2547     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
2548     &                          +zdqcloud(1,:,igcm_h2o_vap))
2549     
2550          call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
2551     &                    rice)
2552         ENDIF ! of IF (water)
2553         
2554ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2555ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2556
2557
2558         zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
2559
2560         do l=2,nlayer-1
2561            tmean=zt(1,l)
2562            if(zt(1,l).ne.zt(1,l-1))
2563     &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
2564              zlocal(l)= zlocal(l-1)
2565     &        -log(zplay(1,l)/zplay(1,l-1))*rnew(1,l)*tmean/g
2566         enddo
2567         zlocal(nlayer)= zlocal(nlayer-1)-
2568     &                   log(zplay(1,nlayer)/zplay(1,nlayer-1))*
2569     &                   rnew(1,nlayer)*tmean/g
2570#endif
2571
2572      END IF       ! if(ngrid.ne.1)
2573
2574      icount=icount+1
2575      RETURN
2576      END
Note: See TracBrowser for help on using the repository browser.