source: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 @ 726

Last change on this file since 726 was 726, checked in by jleconte, 13 years ago

17/07/2012 == JL for LK

  • Generalization of aerosol scheme:
    • any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
    • addition of a module with the id numbers for aerosols (aerosol_mod.F90).
    • initialization of aerosols id numbers in iniaerosol.F90
    • compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
  • For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

  • If starting from an old start file, recreate start file with the q=0 option in newstart.e.
  • update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

  • To add other aerosols, see Laura Kerber.
  • Property svn:executable set to *
File size: 70.1 KB
Line 
1      subroutine physiq(ngrid,nlayer,nq,   &
2                  firstcall,lastcall,      &
3                  pday,ptime,ptimestep,    &
4                  pplev,pplay,pphi,        &
5                  pu,pv,pt,pq,             &
6                  pw,                      &
7                  pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
8 
9      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
10      use watercommon_h, only : RLVTT
11      use gases_h
12      use radcommon_h, only: sigma
13      use aerosol_mod
14      implicit none
15
16
17!==================================================================
18!     
19!     Purpose
20!     -------
21!     Central subroutine for all the physics parameterisations in the
22!     universal model. Originally adapted from the Mars LMDZ model.
23!
24!     The model can be run without or with tracer transport
25!     depending on the value of "tracer" in file "callphys.def".
26!
27!
28!   It includes:
29!
30!      1.  Initialization:
31!      1.1 Firstcall initializations
32!      1.2 Initialization for every call to physiq
33!      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
34!      2. Compute radiative transfer tendencies
35!         (longwave and shortwave).
36!      4. Vertical diffusion (turbulent mixing):
37!      5. Convective adjustment
38!      6. Condensation and sublimation of gases (currently just CO2).
39!      7. TRACERS
40!       7a. water and water ice
41!       7c. other schemes for tracer transport (lifting, sedimentation)
42!       7d. updates (pressure variations, surface budget)
43!      9. Surface and sub-surface temperature calculations
44!     10. Write outputs :
45!           - "startfi", "histfi" if it's time
46!           - Saving statistics if "callstats = .true."
47!           - Output any needed variables in "diagfi"
48!     10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
49!
50!   arguments
51!   ---------
52!
53!   input
54!   -----
55!    ecri                  period (in dynamical timestep) to write output
56!    ngrid                 Size of the horizontal grid.
57!                          All internal loops are performed on that grid.
58!    nlayer                Number of vertical layers.
59!    nq                    Number of advected fields
60!    firstcall             True at the first call
61!    lastcall              True at the last call
62!    pday                  Number of days counted from the North. Spring
63!                          equinoxe.
64!    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
65!    ptimestep             timestep (s)
66!    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
67!    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
68!    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
69!    pu(ngrid,nlayer)      u component of the wind (ms-1)
70!    pv(ngrid,nlayer)      v component of the wind (ms-1)
71!    pt(ngrid,nlayer)      Temperature (K)
72!    pq(ngrid,nlayer,nq)   Advected fields
73!    pudyn(ngrid,nlayer)    \
74!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
75!    ptdyn(ngrid,nlayer)     / corresponding variables
76!    pqdyn(ngrid,nlayer,nq) /
77!    pw(ngrid,?)           vertical velocity
78!
79!   output
80!   ------
81!
82!    pdu(ngrid,nlayermx)        \
83!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
84!    pdt(ngrid,nlayermx)         /  variables due to physical processes.
85!    pdq(ngrid,nlayermx)        /
86!    pdpsrf(ngrid)             /
87!    tracerdyn                 call tracer in dynamical part of GCM ?
88!
89!
90!     Authors
91!     -------
92!           Frederic Hourdin    15/10/93
93!           Francois Forget     1994
94!           Christophe Hourdin  02/1997
95!           Subroutine completely rewritten by F. Forget (01/2000)
96!           Water ice clouds: Franck Montmessin (update 06/2003)
97!           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
98!           New correlated-k radiative scheme: R. Wordsworth (2009)
99!           Many specifically Martian subroutines removed: R. Wordsworth (2009)       
100!           Improved water cycle: R. Wordsworth / B. Charnay (2010)
101!           To F90: R. Wordsworth (2010)
102!           New turbulent diffusion scheme: J. Leconte (2012)
103!           Loops converted to F90 matrix format: J. Leconte (2012)
104!     
105!==================================================================
106
107
108!    0.  Declarations :
109!    ------------------
110
111#include "dimensions.h"
112#include "dimphys.h"
113#include "comgeomfi.h"
114#include "surfdat.h"
115#include "comsoil.h"
116#include "comdiurn.h"
117#include "callkeys.h"
118#include "comcstfi.h"
119#include "planete.h"
120#include "comsaison.h"
121#include "control.h"
122#include "tracer.h"
123#include "watercap.h"
124#include "netcdf.inc"
125
126! Arguments :
127! -----------
128
129!   inputs:
130!   -------
131      integer ngrid,nlayer,nq
132      real ptimestep
133      real pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
134      real pphi(ngridmx,nlayer)
135      real pu(ngridmx,nlayer),pv(ngridmx,nlayer)
136      real pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
137      real pw(ngridmx,nlayer)        ! pvervel transmitted by dyn3d
138      real zh(ngridmx,nlayermx)      ! potential temperature (K)
139      logical firstcall,lastcall
140
141      real pday
142      real ptime
143      logical tracerdyn
144
145!   outputs:
146!   --------
147!     physical tendencies
148      real pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
149      real pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
150      real pdpsrf(ngridmx) ! surface pressure tendency
151
152
153! Local saved variables:
154! ----------------------
155!     aerosol (dust or ice) extinction optical depth  at reference wavelength
156!     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
157!      aerosol optical properties:
158!      real aerosol(ngridmx,nlayermx,naerkind)
159!     this is now internal to callcorrk and hence no longer needed here
160
161      integer day_ini                ! Initial date of the run (sol since Ls=0)
162      integer icount                 ! counter of calls to physiq during the run.
163      real tsurf(ngridmx)            ! Surface temperature (K)
164      real tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
165      real albedo(ngridmx)           ! Surface albedo
166
167      real albedo0(ngridmx)          ! Surface albedo
168      integer rnat(ngridmx)          ! added by BC
169      save rnat
170
171      real emis(ngridmx)             ! Thermal IR surface emissivity
172      real dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
173      real fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
174      real fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
175      real capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
176      real fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
177      real qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
178      real q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
179
180      save day_ini, icount
181      save tsurf,tsoil
182      save albedo0,albedo,emis,q2
183      save capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
184
185! Local variables :
186! -----------------
187
188!     aerosol (dust or ice) extinction optical depth  at reference wavelength
189!     for the "naerkind" optically active aerosols:
190      real aerosol(ngridmx,nlayermx,naerkind)
191
192      character*80 fichier
193      integer l,ig,ierr,iq,i, tapphys,nw
194
195      real fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
196      real fluxsurf_sw(ngridmx)      ! incident SW (stellar) surface flux (W.m-2)
197      real fluxtop_lw(ngridmx)       ! Outgoing LW (IR) flux to space (W.m-2)
198      real fluxabs_sw(ngridmx)       ! Absorbed SW (stellar) flux (W.m-2)
199
200      real fluxtop_dn(ngridmx)
201      real fluxdyn(ngridmx)          ! horizontal heat transport by dynamics     
202      real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
203      real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
204 
205      real,save :: sensibFlux(ngridmx)   ! turbulent flux given by the atm to the surface
206 
207      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu
208
209
210      real zls                       ! solar longitude (rad)
211      real zday                      ! date (time since Ls=0, in martian days)
212      real zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
213      real zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
214      real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
215
216      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
217
218!     Tendencies due to various processes:
219      real dqsurf(ngridmx,nqmx)
220      real,save :: zdtlw(ngridmx,nlayermx)                    ! (K/s)
221      real,save :: zdtsw(ngridmx,nlayermx)                    ! (K/s)
222
223      real cldtlw(ngridmx,nlayermx)                           ! (K/s) LW heating rate for clear areas
224      real cldtsw(ngridmx,nlayermx)                           ! (K/s) SW heating rate for clear areas
225      real zdtsurf(ngridmx)                                   ! (K/s)
226      real dtlscale(ngridmx,nlayermx)
227      real zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
228      real zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
229      real zdtdif(ngridmx,nlayermx)                           ! (K/s)
230      real zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
231      real zdhadj(ngridmx,nlayermx)                           ! (K/s)
232      real zdtgw(ngridmx,nlayermx)                            ! (K/s)
233      real zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
234      real zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
235      real zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
236
237      real zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
238      real zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
239      real zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
240      real zdqadj(ngridmx,nlayermx,nqmx)
241      real zdqc(ngridmx,nlayermx,nqmx)
242      real zdqlscale(ngridmx,nlayermx,nqmx)
243      real zdqslscale(ngridmx,nqmx)
244      real zdqchim(ngridmx,nlayermx,nqmx)
245      real zdqschim(ngridmx,nqmx)
246
247      real zdteuv(ngridmx,nlayermx)    ! (K/s)
248      real zdtconduc(ngridmx,nlayermx) ! (K/s)
249      real zdumolvis(ngridmx,nlayermx)
250      real zdvmolvis(ngridmx,nlayermx)
251      real zdqmoldiff(ngridmx,nlayermx,nqmx)
252
253!     Local variables for local calculations:
254      real zflubid(ngridmx)
255      real zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
256      real zdum1(ngridmx,nlayermx)
257      real zdum2(ngridmx,nlayermx)
258      real ztim1,ztim2,ztim3, z1,z2
259      real ztime_fin
260      real zdh(ngridmx,nlayermx)
261      integer length
262      parameter (length=100)
263
264! local variables only used for diagnostics (output in file "diagfi" or "stats")
265! ------------------------------------------------------------------------------
266      real ps(ngridmx), zt(ngridmx,nlayermx)
267      real zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
268      real zq(ngridmx,nlayermx,nqmx)
269      character*2 str2
270      character*5 str5
271      real zdtadj(ngridmx,nlayermx)
272      real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
273      save ztprevious
274      real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
275      real qtot1,qtot2            ! total aerosol mass
276      integer igmin, lmin
277      logical tdiag
278
279      real zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
280      real zstress(ngrid), cd
281      real hco2(nqmx), tmean, zlocal(nlayermx)
282      real vmr(ngridmx,nlayermx) ! volume mixing ratio
283
284      real time_phys
285
286!     reinstated by RW for diagnostic
287      real tau_col(ngridmx)
288      save tau_col
289     
290!     included by RW to reduce insanity of code
291      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
292
293!     included by RW to compute tracer column densities
294      real qcol(ngridmx,nqmx)
295
296!     included by RW for H2O precipitation
297      real zdtrain(ngridmx,nlayermx)
298      real zdqrain(ngridmx,nlayermx,nqmx)
299      real zdqsrain(ngridmx)
300      real zdqssnow(ngridmx)
301      real rainout(ngridmx)
302
303!     included by RW for H2O Manabe scheme
304      real dtmoist(ngridmx,nlayermx)
305      real dqmoist(ngridmx,nlayermx,nqmx)
306
307      real qvap(ngridmx,nlayermx)
308      real dqvaplscale(ngridmx,nlayermx)
309      real dqcldlscale(ngridmx,nlayermx)
310      real rneb_man(ngridmx,nlayermx)
311      real rneb_lsc(ngridmx,nlayermx)
312
313!     included by RW to account for surface cooling by evaporation
314      real dtsurfh2olat(ngridmx)
315
316
317!     to test energy conservation (RW+JL)
318      real mass(ngridmx,nlayermx),massarea(ngridmx,nlayermx)
319      real dEtot, dEtots, AtmToSurf_TurbFlux
320      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
321      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdiff(ngridmx,nlayermx)
322      real dEdiffs(ngridmx),dEdiff(ngridmx)
323      real madjdE(ngridmx), lscaledE(ngridmx)
324!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
325     
326      real dItot, dVtot
327
328!     included by BC for evaporation     
329      real qevap(ngridmx,nlayermx,nqmx)
330      real tevap(ngridmx,nlayermx)
331      real dqevap(ngridmx,nlayermx)
332      real dtevap(ngridmx,nlayermx)
333
334!     included by BC for hydrology
335      real hice(ngridmx)
336
337!     included by RW to test water conservation (by routine)
338      real h2otot
339      real dWtot, dWtots
340      real h2o_surf_all
341      logical watertest
342      save watertest
343
344!     included by RW for RH diagnostic
345      real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx)
346
347!     included by RW for hydrology
348      real dqs_hyd(ngridmx,nqmx)
349      real zdtsurf_hyd(ngridmx)
350
351!     included by RW for water cycle conservation tests
352      real icesrf,liqsrf,icecol,vapcol
353
354!     included by BC for double radiative transfer call
355      logical clearsky
356      real zdtsw1(ngridmx,nlayermx)
357      real zdtlw1(ngridmx,nlayermx)
358      real fluxsurf_lw1(ngridmx)     
359      real fluxsurf_sw1(ngridmx) 
360      real fluxtop_lw1(ngridmx)   
361      real fluxabs_sw1(ngridmx) 
362      real tau_col1(ngrid)
363      real OLR_nu1(ngrid,L_NSPECTI)
364      real OSR_nu1(ngrid,L_NSPECTV)
365      real tf, ntf
366
367!     included by BC for cloud fraction computations
368      real cloudfrac(ngridmx,nlayermx)
369      real totcloudfrac(ngridmx)
370
371!     included by RW for vdifc water conservation test
372      real nconsMAX
373      real vdifcncons(ngridmx)
374      real cadjncons(ngridmx)
375
376!      double precision qsurf_hist(ngridmx,nqmx)
377      real qsurf_hist(ngridmx,nqmx)
378      save qsurf_hist
379
380!     included by RW for temp convadj conservation test
381      real playtest(nlayermx)
382      real plevtest(nlayermx)
383      real ttest(nlayermx)
384      real qtest(nlayermx)
385      integer igtest
386
387!     included by RW for runway greenhouse 1D study
388      real muvar(ngridmx,nlayermx+1)
389
390!     included by RW for variable H2O particle sizes
391      real reffH2O(ngridmx,nlayermx)
392      real reffcol(ngridmx,naerkind)
393
394!     included by RW for sourceevol
395      real ice_initial(ngridmx)!, isoil
396      real delta_ice,ice_tot
397      real ice_min(ngridmx)
398      save ice_min
399      save ice_initial
400
401      integer num_run
402      logical ice_update
403      save ice_update
404
405!=======================================================================
406     
407
408! 1. Initialisation
409! -----------------
410
411!  1.1   Initialisation only at first call
412!  ---------------------------------------
413      if (firstcall) then
414
415
416!        variables set to 0
417!        ~~~~~~~~~~~~~~~~~~
418         dtrad(:,:) = 0.0
419         fluxrad(:) = 0.0
420         tau_col(:) = 0.0
421         zdtsw(:,:) = 0.0
422         zdtlw(:,:) = 0.0
423
424!        initialize aerosol indexes
425!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
426            call iniaerosol()
427
428     
429!        initialize tracer names, indexes and properties
430!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
431         tracerdyn=tracer
432         if (tracer) then
433            call initracer()
434         endif ! end tracer
435
436!
437
438!        read startfi (initial parameters)
439!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
440         call phyetat0("startfi.nc",0,0,nsoilmx,nq,           &
441               day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
442               cloudfrac,totcloudfrac,hice)
443
444         if (pday.ne.day_ini) then
445           write(*,*) "ERROR in physiq.F90:"
446           write(*,*) "bad synchronization between physics and dynamics"
447           write(*,*) "dynamics day: ",pday
448           write(*,*) "physics day:  ",day_ini
449           stop
450         endif
451
452         write (*,*) 'In physiq day_ini =', day_ini
453
454!        Initialize albedo and orbital calculation
455!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
456         call surfini(ngrid,qsurf,albedo0)
457         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
458
459         do ig=1,ngrid
460            albedo(ig)=albedo0(ig)
461         enddo
462
463         if(tlocked)then
464            print*,'Planet is tidally locked at resonance n=',nres
465            print*,'Make sure you have the right rotation rate!!!'
466         endif
467
468!        initialize soil
469!        ~~~~~~~~~~~~~~~
470         if (callsoil) then
471            call soil(ngrid,nsoilmx,firstcall,inertiedat,     &
472                ptimestep,tsurf,tsoil,capcal,fluxgrd)
473         else
474            print*,'WARNING! Thermal conduction in the soil turned off'
475            do ig=1,ngrid
476               capcal(ig)=1.e6
477               fluxgrd(ig)=0.
478               if(noradsurf)then
479                  fluxgrd(ig)=10.0
480               endif
481            enddo
482            print*,'Flux from ground = ',fluxgrd,' W m^-2'
483         endif
484         icount=1
485
486!        decide whether to update ice at end of run
487!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
488         ice_update=.false.
489         if(sourceevol)then
490            open(128,file='num_run',form='formatted')
491            read(128,*) num_run
492            close(128)
493       
494            if(num_run.ne.0.and.mod(num_run,2).eq.0)then
495            !if(num_run.ne.0.and.mod(num_run,3).eq.0)then
496               print*,'Updating ice at end of this year!'
497               ice_update=.true.
498               ice_min(:)=1.e4
499            endif 
500         endif
501
502!        define surface as continent or ocean
503!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
504         do ig=1,ngridmx
505            rnat(ig)=1
506
507!            if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
508            if(inertiedat(ig,1).gt.1.E4)then
509            rnat(ig)=0
510            endif
511         enddo
512
513         print*,'WARNING! Surface type currently decided by surface inertia'
514         print*,'This should be improved e.g. in newstart.F'
515
516
517!        initialise surface history variable
518!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
519         do ig=1,ngridmx
520            do iq=1,nqmx
521               qsurf_hist(ig,iq)=qsurf(ig,iq)
522            enddo
523         enddo
524
525!        initialise variable for dynamical heating diagnostic
526!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
527         ztprevious(:,:)=pt(:,:)
528
529!        Set temperature just above condensation temperature (for Early Mars)
530!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
531         if(nearco2cond) then
532            write(*,*)' WARNING! Starting at Tcond+1K'
533            do l=1, nlayer
534               do ig=1,ngrid
535                  pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4     &
536                      -pt(ig,l)) / ptimestep
537               enddo
538            enddo
539         endif
540
541         if(meanOLR)then
542            ! to record global radiative balance
543            call system('rm -f rad_bal.out')
544            ! to record global mean/max/min temperatures
545            call system('rm -f tem_bal.out')
546            ! to record global hydrological balance
547            call system('rm -f h2o_bal.out')
548         endif
549
550         watertest=.false.
551         if(water)then
552            ! initialise variables for water cycle
553
554            if(enertest)then
555               watertest = .true.
556            endif
557
558            do ig=1,ngrid
559!already done above (JL12)
560!               qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
561               if(ice_update)then
562                  ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
563               endif
564            enddo
565
566         endif
567         call su_watercycle ! even if we don't have a water cycle, we might
568                            ! need epsi for the wvp definitions in callcorrk.F
569
570      endif        !  (end of "if firstcall")
571
572! ---------------------------------------------------
573! 1.2   Initializations done at every physical timestep:
574! ---------------------------------------------------
575
576      if (ngrid.NE.ngridmx) then
577         print*,'STOP in PHYSIQ'
578         print*,'Probleme de dimensions :'
579         print*,'ngrid     = ',ngrid
580         print*,'ngridmx   = ',ngridmx
581         stop
582      endif
583
584!     Initialize various variables
585!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
586
587      pdu(:,:) = 0.0
588      pdv(:,:) = 0.0
589!      if ( (.not.nearco2cond).and.(.not.firstcall) ) then
590      if ( .not.nearco2cond ) then
591         pdt(:,:) = 0.0
592      endif ! this was the source of an evil bug...
593      pdq(:,:,:)  = 0.0
594      pdpsrf(:)   = 0.0
595      zflubid(:)  = 0.0
596      zdtsurf(:)  = 0.0
597      dqsurf(:,:) = 0.0
598
599      zday=pday+ptime ! compute time, in sols (and fraction thereof)
600
601!     Compute Stellar Longitude (Ls)
602!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
603      if (season) then
604         call stellarlong(zday,zls)
605      else
606         call stellarlong(float(day_ini),zls)
607      end if
608
609!     Compute geopotential between layers
610!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
611
612      do l=1,nlayer
613         do ig=1,ngrid
614            zzlay(ig,l)=pphi(ig,l)/g
615         enddo
616      enddo
617      do ig=1,ngrid
618         zzlev(ig,1)=0.
619         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
620      enddo
621      do l=2,nlayer
622         do ig=1,ngrid
623            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
624            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
625            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
626         enddo
627      enddo
628!     Potential temperature calculation may not be the same in physiq and dynamic...
629
630!     Compute potential temperature
631!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
632
633
634      do l=1,nlayer         
635         do ig=1,ngrid
636            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
637            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
638            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/g
639            massarea(ig,l)=mass(ig,l)*area(ig)
640         enddo
641      enddo
642
643!-----------------------------------------------------------------------
644!    2. Compute radiative tendencies
645!-----------------------------------------------------------------------
646
647      if (callrad) then
648         if( mod(icount-1,iradia).eq.0.or.lastcall) then
649
650!          Compute local stellar zenith angles
651!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
652           call orbite(zls,dist_star,declin)
653
654           if (tlocked) then
655              ztim1=SIN(declin)
656              ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres)
657              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
658
659               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
660               ztim1,ztim2,ztim3,mu0,fract)
661
662           elseif (diurnal) then
663               ztim1=SIN(declin)
664               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
665               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
666
667               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
668               ztim1,ztim2,ztim3,mu0,fract)
669
670           else
671
672               call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
673               ! WARNING: this function appears not to work in 1D
674
675           endif
676
677           if (corrk) then
678
679!          a) Call correlated-k radiative transfer scheme
680!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
681
682              if(kastprof)then
683                 print*,'kastprof should not = true here'
684                 call abort
685              endif
686              muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
687     
688!             standard callcorrk
689              clearsky=.false.
690              call callcorrk(ngrid,nlayer,pq,nq,qsurf,                   &
691                      albedo,emis,mu0,pplev,pplay,pt,                    &
692                      tsurf,fract,dist_star,aerosol,muvar,               &
693                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
694                      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
695                      reffrad,tau_col,cloudfrac,totcloudfrac,            &
696                      clearsky,firstcall,lastcall)     
697
698!             Option to call scheme once more for clear regions
699              if(CLFvarying)then
700
701                 ! ---> PROBLEMS WITH ALLOCATED ARRAYS
702                 ! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
703                 clearsky=.true.
704                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
705                      albedo,emis,mu0,pplev,pplay,pt,                      &
706                      tsurf,fract,dist_star,aerosol,muvar,                 &
707                      zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
708                      fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
709                      reffrad,tau_col1,cloudfrac,totcloudfrac,             &
710                      clearsky,firstcall,lastcall)
711                 clearsky = .false.  ! just in case.     
712
713                 ! Sum the fluxes and heating rates from cloudy/clear cases
714                 do ig=1,ngrid
715                    tf=totcloudfrac(ig)
716                    ntf=1.-tf
717                   
718                    fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
719                    fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
720                    fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
721                    fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
722                    tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
723                   
724                    do l=1,nlayer
725                       zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw(ig,l)
726                       zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw(ig,l)
727                    enddo
728
729                    do nw=1,L_NSPECTV
730                       OSR_nu(ig,nw) = ntf*OSR_nu1(ig,nw) + tf*OSR_nu(ig,nw)                   
731                    enddo
732                    do nw=1,L_NSPECTI
733                       OLR_nu(ig,nw) = ntf*OLR_nu1(ig,nw) + tf*OLR_nu(ig,nw)                   
734                    enddo
735
736                 enddo
737
738              endif !CLFvarying
739
740!             Radiative flux from the sky absorbed by the surface (W.m-2)
741              GSR=0.0
742              do ig=1,ngrid
743                 fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)             &
744                      +fluxsurf_sw(ig)*(1.-albedo(ig))
745
746                 if(noradsurf)then ! no lower surface; SW flux just disappears
747                    GSR = GSR + fluxsurf_sw(ig)*area(ig)
748                    fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
749                 endif
750
751              enddo
752              if(noradsurf)then
753                 print*,'SW lost in deep atmosphere = ',GSR/totarea,' W m^-2'
754              endif
755
756!             Net atmospheric radiative heating rate (K.s-1)
757              do l=1,nlayer
758                 do ig=1,ngrid
759                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
760                 enddo
761              enddo
762
763
764           elseif(newtonian)then
765
766!          b) Call Newtonian cooling scheme
767!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
768              call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
769
770              do ig=1,ngrid
771                 zdtsurf(ig) = +(pt(ig,1)-tsurf(ig))/ptimestep
772              enddo
773              ! e.g. surface becomes proxy for 1st atmospheric layer ?
774
775           else
776
777!          c) Atmosphere has no radiative effect
778!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
779              do ig=1,ngrid
780                 fluxtop_dn(ig)  = fract(ig)*mu0(ig)*Fat1AU/dist_star**2
781                 if(ngrid.eq.1)then ! / by 4 globally in 1D case...
782                    fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
783                 endif
784                 fluxsurf_sw(ig) = fluxtop_dn(ig)
785                 fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
786                 fluxtop_lw(ig)  = emis(ig)*sigma*tsurf(ig)**4
787              enddo ! radiation skips the atmosphere entirely
788
789              do l=1,nlayer
790                 do ig=1,ngrid
791                    dtrad(ig,l)=0.0
792                 enddo
793              enddo ! hence no atmospheric radiative heating
794
795           endif                ! if corrk
796
797        endif ! of if(mod(icount-1,iradia).eq.0)
798
799
800!    Transformation of the radiative tendencies
801!    ------------------------------------------
802
803        do ig=1,ngrid
804           zplanck(ig)=tsurf(ig)*tsurf(ig)
805           zplanck(ig)=emis(ig)*sigma*zplanck(ig)*zplanck(ig)
806           fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
807        enddo
808
809        do l=1,nlayer
810           do ig=1,ngrid
811              pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
812           enddo
813        enddo
814
815!-------------------------
816! test energy conservation
817         if(enertest)then
818            dEtotSW  = cpp*SUM(massarea(:,:)*zdtsw(:,:))/totarea
819            dEtotLW  = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea
820            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea
821            dEtotsLW = SUM((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:))/totarea
822            dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
823            dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
824            print*,'---------------------------------------------------------------'
825            print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
826            print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
827            print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
828            print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
829            print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
830            print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
831         endif
832!-------------------------
833
834      endif ! of if (callrad)
835
836!-----------------------------------------------------------------------
837!    4. Vertical diffusion (turbulent mixing):
838!    -----------------------------------------
839
840      if (calldifv) then
841     
842         do ig=1,ngrid
843            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
844         enddo
845
846         zdum1(:,:)=0.0
847         zdum2(:,:)=0.0
848
849
850!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff     
851         if (UseTurbDiff) then
852         
853           call turbdiff(ngrid,nlayer,nq,rnat,       &
854              ptimestep,capcal,lwrite,               &
855              pplay,pplev,zzlay,zzlev,z0,            &
856              pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
857              zdum1,zdum2,pdt,pdq,zflubid,           &
858              zdudif,zdvdif,zdtdif,zdtsdif,          &
859              sensibFlux,q2,zdqdif,zdqsdif,lastcall)
860
861         else
862         
863           do l=1,nlayer
864             do ig=1,ngrid
865               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
866             enddo
867           enddo
868 
869           call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
870              ptimestep,capcal,lwrite,               &
871              pplay,pplev,zzlay,zzlev,z0,            &
872              pu,pv,zh,pq,tsurf,emis,qsurf,          &
873              zdum1,zdum2,zdh,pdq,zflubid,           &
874              zdudif,zdvdif,zdhdif,zdtsdif,          &
875              sensibFlux,q2,zdqdif,zdqsdif,lastcall)
876
877           do l=1,nlayer
878             do ig=1,ngrid
879                zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
880             enddo
881           enddo
882
883         end if !turbdiff
884
885         do l=1,nlayer
886            do ig=1,ngrid
887               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
888               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
889               pdt(ig,l)=pdt(ig,l)+zdtdif(ig,l)
890            enddo
891         enddo
892
893         do ig=1,ngrid
894            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
895         enddo
896
897         if (tracer) then
898           do iq=1, nq
899            do l=1,nlayer
900              do ig=1,ngrid
901                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
902              enddo
903            enddo
904           enddo
905           do iq=1, nq
906              do ig=1,ngrid
907                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
908              enddo
909           enddo
910
911         end if ! of if (tracer)
912
913         !-------------------------
914         ! test energy conservation
915         if(enertest)then
916            dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
917            do ig = 1, ngrid
918               dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
919               dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
920            enddo
921            dEtot = SUM(dEdiff(:)*area(:))/totarea
922            dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
923            dEtots = SUM(dEdiffs(:)*area(:))/totarea
924            AtmToSurf_TurbFlux=SUM(sensibFlux(:)*area(:))/totarea
925            if (UseTurbDiff) then
926               print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
927               print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
928               print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
929            else
930               print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
931               print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
932               print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
933            end if
934! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
935!    but not given back elsewhere
936         endif
937         !-------------------------
938
939         !-------------------------
940         ! test water conservation
941         if(watertest.and.water)then
942            dWtot = SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap))*ptimestep/totarea
943            dWtots =  SUM(zdqsdif(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
944            do ig = 1, ngrid
945               vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
946            Enddo
947            dWtot = dWtot + SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice))*ptimestep/totarea
948            dWtots = dWtots + SUM(zdqsdif(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
949            do ig = 1, ngrid
950               vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
951            Enddo           
952            nconsMAX=MAXVAL(vdifcncons(:))
953
954            print*,'---------------------------------------------------------------'
955            print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
956            print*,'In difv surface water change            =',dWtots,' kg m-2'
957            print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
958            print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
959
960         endif
961         !-------------------------
962
963      else   
964
965         if(.not.newtonian)then
966
967            do ig=1,ngrid
968               zdtsurf(ig) = zdtsurf(ig) + (fluxrad(ig) + fluxgrd(ig))/capcal(ig)
969            enddo
970
971         endif
972
973      endif ! of if (calldifv)
974
975
976!-----------------------------------------------------------------------
977!   5. Dry convective adjustment:
978!   -----------------------------
979
980      if(calladj) then
981
982         do l=1,nlayer
983            do ig=1,ngrid
984               zdh(ig,l) = pdt(ig,l)/zpopsk(ig,l)
985            enddo
986         enddo
987         zduadj(:,:)=0.0
988         zdvadj(:,:)=0.0
989         zdhadj(:,:)=0.0
990
991
992         call convadj(ngrid,nlayer,nq,ptimestep,    &
993              pplay,pplev,zpopsk,                   &
994              pu,pv,zh,pq,                          &
995              pdu,pdv,zdh,pdq,                      &
996              zduadj,zdvadj,zdhadj,                 &
997              zdqadj)
998
999         do l=1,nlayer
1000            do ig=1,ngrid
1001               pdu(ig,l) = pdu(ig,l) + zduadj(ig,l)
1002               pdv(ig,l) = pdv(ig,l) + zdvadj(ig,l)
1003               pdt(ig,l)    = pdt(ig,l) + zdhadj(ig,l)*zpopsk(ig,l)
1004               zdtadj(ig,l) = zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
1005            enddo
1006         enddo
1007
1008         if(tracer) then
1009           do iq=1, nq
1010            do l=1,nlayer
1011              do ig=1,ngrid
1012                 pdq(ig,l,iq) = pdq(ig,l,iq) + zdqadj(ig,l,iq)
1013              enddo
1014            enddo
1015           enddo
1016         end if
1017
1018         !-------------------------
1019         ! test energy conservation
1020         if(enertest)then
1021            dEtot=cpp*SUM(massarea(:,:)*zdtadj(:,:))/totarea
1022          print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
1023         endif
1024         !-------------------------
1025
1026         !-------------------------
1027         ! test water conservation
1028         if(watertest)then
1029            dWtot = SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap))*ptimestep/totarea
1030            do ig = 1, ngrid
1031               cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
1032            Enddo
1033            dWtot = dWtot + SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice))*ptimestep/totarea
1034            do ig = 1, ngrid
1035               cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
1036            Enddo           
1037            nconsMAX=MAXVAL(cadjncons(:))
1038
1039            print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
1040            print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
1041         endif
1042         !-------------------------
1043               
1044      endif ! of if(calladj)
1045
1046!-----------------------------------------------------------------------
1047!   6. Carbon dioxide condensation-sublimation:
1048!   -------------------------------------------
1049
1050      if (co2cond) then
1051         if (.not.tracer) then
1052            print*,'We need a CO2 ice tracer to condense CO2'
1053            call abort
1054         endif
1055     print*, 'we are in co2cond!!!!!'
1056         call condense_cloud(ngrid,nlayer,nq,ptimestep,   &
1057              capcal,pplay,pplev,tsurf,pt,                  &
1058              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
1059              qsurf(1,igcm_co2_ice),albedo,emis,            &
1060              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
1061              zdqc,reffrad)
1062
1063         do l=1,nlayer
1064           do ig=1,ngrid
1065             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
1066             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
1067             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
1068           enddo
1069         enddo
1070         do ig=1,ngrid
1071            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
1072         enddo
1073
1074         do iq=1,nq ! should use new notation here !
1075            do l=1,nlayer
1076               do ig=1,ngrid
1077                  pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
1078               enddo
1079            enddo
1080         enddo
1081         ! Note: we do not add surface co2ice tendency
1082         ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
1083
1084         !-------------------------
1085         ! test energy conservation
1086         if(enertest)then
1087            dEtot = cpp*SUM(massarea(:,:)*zdtc(:,:))/totarea
1088            dEtots = SUM(capcal(:)*zdtsurfc(:)*area(:))/totarea
1089            print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
1090            print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
1091         endif
1092         !-------------------------
1093
1094      endif  ! of if (co2cond)
1095
1096
1097!-----------------------------------------------------------------------
1098!   7. Specific parameterizations for tracers
1099!   -----------------------------------------
1100
1101      if (tracer) then
1102
1103!   7a. Water and ice
1104!     ---------------
1105         if (water) then
1106
1107!        ----------------------------------------
1108!        Water ice condensation in the atmosphere
1109!        ----------------------------------------
1110            if(watercond)then
1111
1112               if(RLVTT.gt.1.e-8)then
1113
1114               ! Re-evaporate cloud water/ice
1115               call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
1116               DO l = 1, nlayer 
1117                  DO ig = 1, ngrid
1118                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
1119                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
1120                     pdt(ig,l)              = pdt(ig,l)+dtevap(ig,l)
1121                  enddo
1122               enddo
1123
1124               call moistadj(pt,qevap,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
1125               do l=1,nlayer
1126                  do ig=1,ngrid
1127                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqmoist(ig,l,igcm_h2o_vap)
1128                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqmoist(ig,l,igcm_h2o_ice)
1129                     pdt(ig,l)              = pdt(ig,l)+dtmoist(ig,l)
1130                  enddo
1131               enddo
1132
1133               !-------------------------
1134               ! test energy conservation
1135               if(enertest)then
1136                  dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
1137                  do ig = 1, ngrid
1138                     madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
1139                  enddo
1140                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
1141                 
1142                ! test energy conservation
1143                  dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
1144                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
1145                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
1146               endif
1147               !-------------------------
1148               
1149
1150               endif
1151               
1152
1153               ! Re-evaporate cloud water/ice
1154               call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
1155               do l = 1, nlayer 
1156                  do ig = 1, ngrid
1157                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
1158                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
1159                     pdt(ig,l)              = pdt(ig,l)+dtevap(ig,l)
1160                  enddo
1161               enddo ! note: we use qevap but not tevap in largescale/moistadj
1162                     ! otherwise is a big mess
1163
1164               call largescale(ptimestep,pplev,pplay,pt,qevap,        & ! a bug was here!
1165                    pdt,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc,reffH2O)
1166               do l=1,nlayer
1167                  do ig=1,ngrid
1168                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqvaplscale(ig,l)
1169                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l)
1170                     pdt(ig,l)              = pdt(ig,l)+dtlscale(ig,l)
1171                     
1172                     if(aeroh2o.and.(.not.aerofixh2o)) then
1173                        reffrad(ig,l,iaero_h2o)=reffH2O(ig,l)
1174                     endif
1175
1176                  enddo
1177               enddo
1178
1179               !-------------------------
1180               ! test energy conservation
1181               if(enertest)then
1182                  dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
1183                  do ig = 1, ngrid
1184                     madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
1185                  enddo
1186                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
1187                 
1188                ! test energy conservation
1189                  dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
1190                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
1191                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
1192               endif
1193               !-------------------------
1194
1195               ! compute cloud fraction
1196               do l = 1, nlayer
1197                  do ig = 1,ngrid
1198                     cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
1199                  enddo
1200               enddo
1201
1202               ! compute total cloud fraction in column
1203               call totalcloudfrac(cloudfrac,totcloudfrac)
1204
1205           endif                ! of if (watercondense)
1206           
1207
1208!        --------------------------------
1209!        Water ice / liquid precipitation
1210!        --------------------------------
1211           if(waterrain)then
1212
1213              zdqrain(:,:,:) = 0.0
1214              zdqsrain(:)    = 0.0
1215              zdqssnow(:)    = 0.0
1216
1217              call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
1218                   zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
1219
1220              do l=1,nlayer
1221                 do ig=1,ngrid
1222                    pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+zdqrain(ig,l,igcm_h2o_vap)
1223                    pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+zdqrain(ig,l,igcm_h2o_ice)
1224                    pdt(ig,l)              = pdt(ig,l)+zdtrain(ig,l)
1225                 enddo
1226              enddo
1227
1228              do ig=1,ngrid
1229                 dqsurf(ig,igcm_h2o_vap) = dqsurf(ig,igcm_h2o_vap)+zdqsrain(ig) ! a bug was here
1230                 dqsurf(ig,igcm_h2o_ice) = dqsurf(ig,igcm_h2o_ice)+zdqssnow(ig)
1231                 rainout(ig)             = zdqsrain(ig)+zdqssnow(ig) ! diagnostic   
1232              enddo
1233
1234
1235
1236               !-------------------------
1237               ! test energy conservation
1238               if(enertest)then
1239                  dEtot=cpp*SUM(massarea(:,:)*zdtrain(:,:))/totarea
1240                 print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
1241                  dItot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))/totarea*RLVTT/cpp
1242                  dItot = dItot + SUM(area(:)*zdqssnow(:))/totarea*RLVTT/cpp
1243                  dVtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
1244                  dVtot = dItot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
1245                  dEtot = dItot + dVtot
1246                 print*,'In rain dItot =',dItot,' W m-2'
1247                 print*,'In rain dVtot =',dVtot,' W m-2'
1248                 print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
1249
1250               ! test water conservation
1251                  dWtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
1252                  dWtot = dWtot + SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))*ptimestep/totarea
1253                  dWtots =  SUM((zdqsrain(:)+zdqssnow(:))*area(:))*ptimestep/totarea
1254                 print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
1255                 print*,'In rain surface water change            =',dWtots,' kg m-2'
1256                 print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
1257              endif
1258              !-------------------------
1259
1260           end if                 ! of if (waterrain)
1261        end if                    ! of if (water)
1262
1263
1264!   7c. Aerosol particles
1265!     -------------------
1266!        -------------
1267!        Sedimentation
1268!        -------------
1269        if (sedimentation) then
1270           zdqsed(:,:,:) = 0.0
1271           zdqssed(:,:)  = 0.0
1272
1273
1274           !-------------------------
1275           ! find qtot
1276           if(watertest)then
1277              iq=3
1278              dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
1279              dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
1280              print*,'Before sedim pq  =',dWtot,' kg m-2'
1281              print*,'Before sedim pdq =',dWtots,' kg m-2'
1282           endif
1283           !-------------------------
1284
1285           call callsedim(ngrid,nlayer,ptimestep,           &
1286                pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
1287
1288           !-------------------------
1289           ! find qtot
1290           if(watertest)then
1291              iq=3
1292              dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
1293              dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
1294              print*,'After sedim pq  =',dWtot,' kg m-2'
1295              print*,'After sedim pdq =',dWtots,' kg m-2'
1296           endif
1297           !-------------------------
1298
1299           do iq=1,nq
1300              ! for now, we only allow H2O ice to sediment
1301              ! and as in rain.F90, whether it falls as rain or snow depends
1302              ! only on the surface temperature
1303              do ig=1,ngrid
1304                 do l=1,nlayer
1305                    pdq(ig,l,iq) = pdq(ig,l,iq) + zdqsed(ig,l,iq)
1306                 enddo
1307                 dqsurf(ig,iq) = dqsurf(ig,iq) + zdqssed(ig,iq)
1308              enddo
1309           enddo
1310
1311           !-------------------------
1312           ! test water conservation
1313           if(watertest)then
1314              dWtot = SUM(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice)))*ptimestep/totarea
1315              dWtots = SUM((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:))*ptimestep/totarea
1316              print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
1317              print*,'In sedim surface ice change             =',dWtots,' kg m-2'
1318              print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
1319           endif
1320           !-------------------------
1321
1322        end if   ! of if (sedimentation)
1323
1324
1325!   7d. Updates
1326!     ---------
1327
1328!       ---------------------------------
1329!       Updating tracer budget on surface
1330!       ---------------------------------
1331
1332         if(hydrology)then
1333
1334            call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,    &
1335               capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
1336            ! note: for now, also changes albedo in the subroutine
1337
1338            do ig=1,ngrid
1339               zdtsurf(ig) = zdtsurf(ig) + zdtsurf_hyd(ig)
1340               do iq=1,nq
1341                  qsurf(ig,iq) = qsurf(ig,iq)+ptimestep*dqs_hyd(ig,iq)
1342               enddo
1343            enddo
1344            ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
1345
1346            !-------------------------
1347            ! test energy conservation
1348            if(enertest)then
1349               dEtots = SUM(area(:)*capcal(:)*zdtsurf_hyd(:))/totarea
1350               print*,'In hydrol surface energy change     =',dEtots,' W m-2'
1351            endif
1352            !-------------------------
1353
1354            !-------------------------
1355            ! test water conservation
1356            if(watertest)then
1357               dWtots =  SUM(dqs_hyd(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
1358               print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
1359               dWtots =  SUM(dqs_hyd(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
1360               print*,'In hydrol surface water change          =',dWtots,' kg m-2'
1361               print*,'---------------------------------------------------------------'
1362            endif
1363            !-------------------------
1364
1365         ELSE     ! of if (hydrology)
1366
1367            do iq=1,nq
1368               do ig=1,ngrid
1369                  qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
1370               enddo
1371            enddo           
1372
1373         END IF   ! of if (hydrology)
1374
1375         ! Add qsurf to qsurf_hist, which is what we save in
1376         ! diagfi.nc etc. At the same time, we set the water
1377         ! content of ocean gridpoints back to zero, in order
1378         ! to avoid rounding errors in vdifc, rain
1379         qsurf_hist(:,:) = qsurf(:,:)
1380
1381         if(ice_update)then
1382            do ig = 1, ngrid
1383               ice_min(ig)=min(ice_min(ig),qsurf(ig,igcm_h2o_ice))
1384            enddo
1385         endif
1386
1387      endif   !  of if (tracer)
1388
1389!-----------------------------------------------------------------------
1390!   9. Surface and sub-surface soil temperature
1391!-----------------------------------------------------------------------
1392
1393
1394!     Increment surface temperature
1395      do ig=1,ngrid
1396         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
1397      enddo
1398
1399!     Compute soil temperatures and subsurface heat flux
1400      if (callsoil) then
1401         call soil(ngrid,nsoilmx,.false.,inertiedat,        &
1402                ptimestep,tsurf,tsoil,capcal,fluxgrd)
1403      endif
1404
1405!-------------------------
1406! test energy conservation
1407      if(enertest)then
1408         dEtots=0.0
1409         do ig = 1, ngrid
1410            dEtots = dEtots + capcal(ig)*zdtsurf(ig)*area(ig)
1411         enddo
1412         dEtots=dEtots/totarea
1413         print*,'Surface energy change                 =',dEtots,' W m-2'
1414      endif
1415!-------------------------
1416
1417!-----------------------------------------------------------------------
1418!  10. Perform diagnostics and write output files
1419!-----------------------------------------------------------------------
1420
1421!    -------------------------------
1422!    Dynamical fields incrementation
1423!    -------------------------------
1424!    For output only: the actual model integration is performed in the dynamics
1425 
1426      ! temperature, zonal and meridional wind
1427      do l=1,nlayer
1428        do ig=1,ngrid
1429          zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
1430          zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
1431          zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
1432
1433          ! diagnostic
1434          zdtdyn(ig,l)     = ztprevious(ig,l)-pt(ig,l)
1435          ztprevious(ig,l) = zt(ig,l)
1436        enddo
1437      enddo
1438
1439      if(firstcall)then
1440         zdtdyn(:,:)=0.0
1441      endif
1442
1443      ! dynamical heating diagnostic
1444      fluxdyn(:)=0.
1445      do ig=1,ngrid
1446         do l=1,nlayer
1447            fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)   &
1448                 *(pplev(ig,l)-pplev(ig,l+1))*cpp/g
1449         enddo
1450      enddo
1451
1452      ! tracers
1453      do iq=1, nq
1454        do l=1,nlayer
1455          do ig=1,ngrid
1456            zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
1457          enddo
1458        enddo
1459      enddo
1460
1461      ! surface pressure
1462      do ig=1,ngrid
1463          ps(ig) = pplev(ig,1) + pdpsrf(ig)*ptimestep
1464      enddo
1465
1466      ! pressure
1467      do l=1,nlayer
1468        do ig=1,ngrid
1469             zplev(ig,l) = pplev(ig,l)/pplev(ig,1)*ps(ig)
1470             zplay(ig,l) = pplay(ig,l)/pplev(ig,1)*ps(ig)
1471        enddo
1472      enddo
1473
1474!     ---------------------------------------------------------
1475!     Surface and soil temperature information
1476!     ---------------------------------------------------------
1477
1478      Ts1 = SUM(area(:)*tsurf(:))/totarea
1479      Ts2 = MINVAL(tsurf(:))
1480      Ts3 = MAXVAL(tsurf(:))
1481      if(callsoil)then
1482         TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
1483         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
1484         print*,Ts1,Ts2,Ts3,TsS
1485      else
1486         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
1487         print*,Ts1,Ts2,Ts3
1488      endif
1489
1490!     ---------------------------------------------------------
1491!     Check the energy balance of the simulation during the run
1492!     ---------------------------------------------------------
1493
1494      if(corrk)then
1495
1496         ISR = SUM(area(:)*fluxtop_dn(:))/totarea
1497         ASR = SUM(area(:)*fluxabs_sw(:))/totarea
1498         OLR = SUM(area(:)*fluxtop_lw(:))/totarea
1499         GND = SUM(area(:)*fluxgrd(:))/totarea
1500         DYN = SUM(area(:)*fluxdyn(:))/totarea
1501         do ig=1,ngrid             
1502            if(fluxtop_dn(ig).lt.0.0)then
1503               print*,'fluxtop_dn has gone crazy'
1504               print*,'fluxtop_dn=',fluxtop_dn(ig)
1505               print*,'tau_col=',tau_col(ig)
1506               print*,'aerosol=',aerosol(ig,:,:)
1507               print*,'temp=   ',pt(ig,:)
1508               print*,'pplay=  ',pplay(ig,:)
1509               call abort
1510            endif
1511         end do
1512                     
1513         if(ngridmx.eq.1)then
1514            DYN=0.0
1515         endif
1516
1517         print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
1518         print*, ISR,ASR,OLR,GND,DYN
1519         
1520         if(enertest)then
1521            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
1522            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
1523            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
1524         endif
1525
1526         if(meanOLR)then
1527            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
1528               ! to record global radiative balance
1529               open(92,file="rad_bal.out",form='formatted',position='append')
1530               write(92,*) zday,ISR,ASR,OLR
1531               close(92)
1532               open(93,file="tem_bal.out",form='formatted',position='append')
1533               write(93,*) zday,Ts1,Ts2,Ts3,TsS
1534               close(93)
1535            endif
1536         endif
1537
1538      endif
1539
1540
1541!     ------------------------------------------------------------------
1542!     Diagnostic to test radiative-convective timescales in code
1543!     ------------------------------------------------------------------
1544      if(testradtimes)then
1545         open(38,file="tau_phys.out",form='formatted',position='append')
1546         ig=1
1547         do l=1,nlayer
1548            write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
1549         enddo
1550         close(38)
1551         print*,'As testradtimes enabled,'
1552         print*,'exiting physics on first call'
1553         call abort
1554      endif
1555
1556!     ---------------------------------------------------------
1557!     Compute column amounts (kg m-2) if tracers are enabled
1558!     ---------------------------------------------------------
1559      if(tracer)then
1560         qcol(:,:)=0.0
1561         do iq=1,nq
1562            do ig=1,ngrid
1563               do l=1,nlayer
1564                  qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *    &
1565                            (pplev(ig,l) - pplev(ig,l+1)) / g
1566               enddo
1567            enddo
1568         enddo
1569
1570         ! Generalised for arbitrary aerosols now. (LK)
1571         reffcol(:,:)=0.0
1572         do ig=1,ngrid
1573            do l=1,nlayer
1574               if(co2cond.and.(iaero_co2.ne.0)) then
1575                  reffcol(ig,iaero_co2) = reffcol(ig,iaero_co2) + zq(ig,l,igcm_co2_ice) * &
1576                                  reffrad(ig,l,iaero_co2) *    &
1577                                  (pplev(ig,l) - pplev(ig,l+1)) / g
1578               endif
1579               if(water.and.(iaero_h2o.ne.0)) then
1580                  reffcol(ig,iaero_h2o) = reffcol(ig,iaero_h2o) + zq(ig,l,igcm_h2o_ice) * &
1581                                  reffrad(ig,l,iaero_h2o) *    &
1582                                  (pplev(ig,l) - pplev(ig,l+1)) / g
1583               endif
1584            enddo
1585         enddo
1586
1587      endif
1588
1589!     ---------------------------------------------------------
1590!     Test for water conservation if water is enabled
1591!     ---------------------------------------------------------
1592
1593      if(water)then
1594
1595         icesrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_ice))/totarea
1596         liqsrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_vap))/totarea
1597         icecol = SUM(area(:)*qcol(:,igcm_h2o_ice))/totarea
1598         vapcol = SUM(area(:)*qcol(:,igcm_h2o_vap))/totarea
1599
1600         h2otot = icesrf + liqsrf + icecol + vapcol
1601
1602         print*,' Total water amount [kg m^-2]: ',h2otot
1603         print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
1604         print*, icesrf,liqsrf,icecol,vapcol
1605
1606         if(meanOLR)then
1607            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
1608               ! to record global water balance
1609               open(98,file="h2o_bal.out",form='formatted',position='append')
1610               write(98,*) zday,icesrf,liqsrf,icecol,vapcol
1611               close(98)
1612            endif
1613         endif
1614
1615      endif
1616
1617!     ---------------------------------------------------------
1618!     Calculate RH for diagnostic if water is enabled
1619!     ---------------------------------------------------------
1620
1621      if(water)then
1622         do l = 1, nlayer
1623            do ig = 1, ngrid
1624               call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
1625               RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
1626            enddo
1627         enddo
1628
1629         ! compute maximum possible H2O column amount (100% saturation)
1630         do ig=1,ngrid
1631            H2Omaxcol(ig)=0.0
1632            do l=1,nlayer
1633               H2Omaxcol(ig) = H2Omaxcol(ig) + qsat(ig,l) *    &
1634                    (pplev(ig,l) - pplev(ig,l+1))/g
1635            enddo
1636         enddo
1637
1638      endif
1639
1640
1641         print*,''
1642         print*,'--> Ls =',zls*180./pi
1643!        -------------------------------------------------------------------
1644!        Writing NetCDF file  "RESTARTFI" at the end of the run
1645!        -------------------------------------------------------------------
1646!        Note: 'restartfi' is stored just before dynamics are stored
1647!              in 'restart'. Between now and the writting of 'restart',
1648!              there will have been the itau=itau+1 instruction and
1649!              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
1650!              thus we store for time=time+dtvr
1651
1652         if(lastcall) then
1653            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
1654
1655
1656!           Update surface ice distribution to iterate to steady state if requested
1657            if(ice_update)then
1658
1659               do ig = 1, ngrid
1660
1661                  delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
1662
1663                  ! add multiple years of evolution
1664                  qsurf_hist(ig,igcm_h2o_ice) = &
1665                     !qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0
1666                     qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
1667
1668                  ! if ice has gone -ve, set to zero
1669                  if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
1670                     qsurf_hist(ig,igcm_h2o_ice) = 0.0
1671                     !qsurf_hist(ig,igcm_h2o_vap) = 0.0
1672                  endif
1673
1674                  ! if ice is seasonal, set to zero (NEW)
1675                  if(ice_min(ig).lt.0.01)then
1676                     qsurf_hist(ig,igcm_h2o_ice) = 0.0
1677                     !qsurf_hist(ig,igcm_h2o_vap) = 0.0
1678                  endif
1679
1680               enddo
1681
1682               ! enforce ice conservation
1683               ice_tot=0.0
1684               do ig = 1, ngrid
1685                  ice_tot = ice_tot + qsurf_hist(ig,igcm_h2o_ice)*area(ig)
1686               enddo
1687               do ig = 1, ngrid
1688                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice)*(icesrf/ice_tot)
1689               enddo
1690
1691            endif
1692
1693            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
1694            call physdem1("restartfi.nc",long,lati,nsoilmx,nq,            &
1695                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
1696                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
1697                    cloudfrac,totcloudfrac,hice)
1698         endif
1699
1700!        -----------------------------------------------------------------
1701!        Saving statistics :
1702!        -----------------------------------------------------------------
1703!        ("stats" stores and accumulates 8 key variables in file "stats.nc"
1704!        which can later be used to make the statistic files of the run:
1705!        "stats")          only possible in 3D runs !
1706
1707         
1708         if (callstats) then
1709
1710            call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
1711            call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1712            call wstats(ngrid,"fluxsurf_lw",                               &
1713                        "Thermal IR radiative flux to surface","W.m-2",2,  &
1714                         fluxsurf_lw)
1715!            call wstats(ngrid,"fluxsurf_sw",                               &
1716!                        "Solar radiative flux to surface","W.m-2",2,       &
1717!                         fluxsurf_sw_tot)
1718            call wstats(ngrid,"fluxtop_lw",                                &
1719                        "Thermal IR radiative flux to space","W.m-2",2,    &
1720                        fluxtop_lw)
1721!            call wstats(ngrid,"fluxtop_sw",                                &
1722!                        "Solar radiative flux to space","W.m-2",2,         &
1723!                        fluxtop_sw_tot)
1724
1725            call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1726            call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1727            call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1728
1729            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
1730            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
1731            call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
1732            call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
1733            call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
1734
1735           if (tracer) then
1736             do iq=1,nq
1737                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1738                call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
1739                          'kg m^-2',2,qsurf(1,iq) )
1740
1741                call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
1742                          'kg m^-2',2,qcol(1,iq) )
1743!                call wstats(ngridmx,trim(noms(iq))//'_reff',                          &
1744!                          trim(noms(iq))//'_reff',                                   &
1745!                          'm',3,reffrad(1,1,iq))
1746              end do
1747            if (water) then
1748               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
1749               call wstats(ngrid,"vmr_h2ovapor",                           &
1750                          "H2O vapour volume mixing ratio","mol/mol",       &
1751                          3,vmr)
1752            endif ! of if (water)
1753
1754           endif !tracer
1755
1756            if(lastcall) then
1757              write (*,*) "Writing stats..."
1758              call mkstats(ierr)
1759            endif
1760          endif !if callstats
1761
1762
1763!        ----------------------------------------------------------------------
1764!        output in netcdf file "DIAGFI", containing any variable for diagnostic
1765!        (output with  period "ecritphy", set in "run.def")
1766!        ----------------------------------------------------------------------
1767!        writediagfi can also be called from any other subroutine for any variable.
1768!        but its preferable to keep all the calls in one place...
1769
1770        call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
1771        call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
1772        call writediagfi(ngrid,"temp","temperature","K",3,zt)
1773        call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
1774        call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
1775        call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
1776        call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
1777        call writediagfi(ngrid,'p','Pressure','Pa',3,pplay)
1778
1779!     Total energy balance diagnostics
1780        if(callrad.and.(.not.newtonian))then
1781           call writediagfi(ngrid,'ALB','Surface albedo',' ',2,albedo)
1782           call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
1783           call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
1784           call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
1785           call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
1786           call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
1787        endif
1788       
1789        if(enertest) then
1790          if (calldifv) then
1791             call writediagfi(ngridmx,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
1792             call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
1793             call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
1794             call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
1795             call writediagfi(ngridmx,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
1796          endif
1797          if (corrk) then
1798             call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
1799             call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
1800          endif
1801          if(watercond) then
1802             call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
1803             call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
1804             call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
1805             call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
1806          endif
1807        endif
1808
1809!     Temporary inclusions for heating diagnostics
1810!        call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
1811!        call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
1812!        call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
1813!        call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
1814
1815        ! debugging
1816        !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
1817        !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
1818
1819!     Output aerosols
1820        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0)   &
1821           call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3, &
1822           reffrad(1,1,iaero_co2))
1823        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0)   &
1824           call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3, &
1825           reffrad(1,1,iaero_h2o))
1826        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0)   &
1827           call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol', &
1828           'um kg m^-2',2,reffcol(1,iaero_co2))
1829        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0)   &
1830           call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol', &
1831           'um kg m^-2',2,reffcol(1,iaero_h2o))
1832
1833!     Output tracers
1834       if (tracer) then
1835        do iq=1,nq
1836          call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
1837!          call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
1838!                          'kg m^-2',2,qsurf(1,iq) )
1839          call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
1840                          'kg m^-2',2,qsurf_hist(1,iq) )
1841          call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
1842                          'kg m^-2',2,qcol(1,iq) )
1843
1844          if(watercond.or.CLFvarying)then
1845!             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
1846!             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
1847!             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
1848             call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
1849          endif
1850
1851          if(waterrain)then
1852             call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
1853             call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
1854          endif
1855
1856          if(hydrology)then
1857             call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
1858          endif
1859
1860          if(ice_update)then
1861             call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
1862             call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
1863          endif
1864
1865          do ig=1,ngrid
1866             if(tau_col(ig).gt.1.e3)then
1867                print*,'WARNING: tau_col=',tau_col(ig)
1868                print*,'at ig=',ig,'in PHYSIQ'
1869             endif
1870          end do
1871
1872          call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
1873
1874        enddo
1875       endif
1876
1877!      output spectrum
1878      if(specOLR.and.corrk)then     
1879                 call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
1880                 call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
1881      endif
1882
1883
1884      icount=icount+1
1885
1886      ! deallocate gas variables
1887      if (lastcall) then
1888        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
1889        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  ! both allocated in su_gases.F90
1890      endif
1891
1892
1893      return
1894    end subroutine physiq
Note: See TracBrowser for help on using the repository browser.