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

Last change on this file since 544 was 544, checked in by acolaitis, 13 years ago

27/02/12 == AC

Continuation of thermals setting, comparisons with mesoscale results on Case C
Added possibility to call gcm (or 1d) with constant prescribed sensible heat flux, in the spirit of direct comparisons with LES

... This is directly comparable to the variable tke_heat_flux in namelist.input
... Requires the use of yamada4.F from terrestrial GCM (mixes more, seem more numerically stable)
... Usually requires high timesteps (>1000) to avoids crashes. Best approach is to compare
height of first model level z1 and teta_1 between LES and 1D, and increase the timesteps until results
between the two models are comparable (might require a slitghly different tke_heat_flux between the two models
due to difference in vertical diffusion schemes and subgrid effects)
... Syntax for use is to add "tke_heat_flux = ..." in callphys.def

Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...)

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