source: trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F @ 1884

Last change on this file since 1884 was 1880, checked in by mvals, 8 years ago

Mars GCM
-In the sub-grid scale cloud scheme zt is replaced is by ztclf and zq by zqclf ('clf' for 'cloud fraction') in order to avoid any confusion for the further schemes, which need initialization.
-zteff and pqeff are initialized in the first part of the CLFvarying scheme, section 0 of the code, instead of been initialized in section 1 (tendencies) with the sub-timesteps.
-The cloud fraction cannot be lower than the settled value mincloud.(MV)

File size: 22.8 KB
RevLine 
[1711]1       MODULE watercloud_mod
2
3       IMPLICIT NONE
4
5       CONTAINS
6
[633]7       SUBROUTINE watercloud(ngrid,nlay,ptimestep,
8     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
[626]9     &                pq,pdq,pdqcloud,pdtcloud,
[358]10     &                nq,tau,tauscaling,rdust,rice,nuice,
[1711]11     &                rsedcloud,rhocloud,totcloudfrac)
[633]12! to use  'getin'
13      USE ioipsl_getincom
[740]14      USE updaterad
[1226]15      USE comcstfi_h
[1036]16      use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
17     &                      igcm_dust_mass, igcm_dust_number,
18     &                      igcm_ccn_mass, igcm_ccn_number,
19     &                      rho_dust, nuice_sed, nuice_ref
[1246]20      use dimradmars_mod, only: naerkind
[38]21      IMPLICIT NONE
22
[633]23
[38]24c=======================================================================
[358]25c  Water-ice cloud formation
26
27c  Includes two different schemes:
28c    - A simplified scheme (see simpleclouds.F)
29c    - An improved microphysical scheme (see improvedclouds.F)
[38]30c
[633]31c  There is a time loop specific to cloud formation
32c  due to timescales smaller than the GCM integration timestep.
33c
[358]34c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
[522]35c           J.-B. Madeleine, Thomas Navarro
[38]36c
[633]37c  2004 - 2012
[38]38c=======================================================================
39
40c-----------------------------------------------------------------------
41c   declarations:
42c   -------------
43
44#include "callkeys.h"
45
46c   Inputs:
47c   ------
48
49      INTEGER ngrid,nlay
[633]50      INTEGER nq                 ! nombre de traceurs
[38]51      REAL ptimestep             ! pas de temps physique (s)
52      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
53      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
[633]54      REAL pdpsrf(ngrid)         ! tendence surf pressure
[38]55      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
56      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
[633]57      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
[38]58
59      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
[633]60      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
[38]61
[1047]62      REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
63      REAL tauscaling(ngrid)   ! Convertion factor for dust amount
64      real rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
[38]65
66c   Outputs:
67c   -------
68
[633]69      real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
70      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
71                                   ! a la chaleur latente
[38]72
73      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
74                               ! (r_c in montmessin_2004)
75      REAL nuice(ngrid,nlay)   ! Estimated effective variance
76                               !   of the size distribution
[1047]77      real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
78      real rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
[38]79
[1711]80      REAL, INTENT(INOUT):: totcloudfrac(ngrid) ! Cloud fraction (A. Pottier 2013)
[38]81c   local:
82c   ------
[633]83     
84      ! for ice radius computation
85      REAL Mo,No
86      REAl ccntyp
87     
88      ! for time loop
89      INTEGER microstep  ! time subsampling step variable
90      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
91      SAVE imicro
92      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
93      SAVE microtimestep
94     
95      ! tendency given by clouds (inside the micro loop)
96      REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
97      REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
[38]98
[633]99      ! global tendency (clouds+physics)
100      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
101      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
102
[1467]103      ! no supersaturation when option supersat is false
104      REAL zt(ngrid,nlay)       ! local value of temperature
105      REAL zqsat(ngrid,nlay)    ! saturation
[633]106
107      INTEGER iq,ig,l
[38]108      LOGICAL,SAVE :: firstcall=.true.
109
[1711]110! Representation of sub-grid water ice clouds A. Pottier 2013
[1880]111      REAL :: ztclf(ngrid, nlay)
112      REAL :: zqclf(ngrid, nlay,nq)
[1711]113      REAL :: zdelt 
114      REAL :: norm
115      REAL :: ponder
116      REAL :: tcond(ngrid,nlay)
[1880]117      REAL :: zqvap(ngrid,nlay)
[1711]118      REAL :: zqice(ngrid,nlay)
119      REAL ::  spant ! delta T for the temperature distribution
120!      REAL :: zqsat(ngrid,nlay) ! saturation
121      REAL :: zteff(ngrid, nlay)! effective temperature in the cloud,neb
122      REAL :: pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud
123      REAL :: cloudfrac(ngrid,nlay) ! cloud fraction
124      REAL :: mincloud ! min cloud frac
125      LOGICAL, save :: flagcloud=.true.
[38]126c    ** un petit test de coherence
127c       --------------------------
128
129      IF (firstcall) THEN
130         
131        if (nq.gt.nqmx) then
132           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
133           write(*,*) 'nq=',nq,' nqmx=',nqmx
134           stop
135        endif
136         
[358]137        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
138        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
[633]139               
140        write(*,*) "time subsampling for microphysic ?"
141#ifdef MESOSCALE
142        imicro = 2
143#else
[951]144        imicro = 30
[633]145#endif
146        call getin("imicro",imicro)
147        write(*,*)"imicro = ",imicro
148       
[38]149        firstcall=.false.
150      ENDIF ! of IF (firstcall)
[1774]151
152      !! AS: moved out of firstcall to allow nesting+evoluting timestep
153      !!     TBD: consider possible diff imicro with domains?
154      microtimestep = ptimestep/real(imicro)
155      write(*,*)"Physical timestep is",ptimestep
156      write(*,*)"Microphysics timestep is",microtimestep
157
[522]158     
[633]159c-----Initialization
160      subpdq(1:ngrid,1:nlay,1:nq) = 0
161      subpdt(1:ngrid,1:nlay)      = 0
162     
163      ! default value if no ice
164      rhocloud(1:ngrid,1:nlay) = rho_dust
[38]165
[1711]166c-------------------------------------------------------------------
167c   0.  Representation of sub-grid water ice clouds
168c------------------
[1880]169c-----Initialization
170      zteff(1:ngrid,1:nlay) = 0
171      pqeff(1:ngrid,1:nlay,1:nq) = 0
172      DO l=1,nlay
173        DO ig=1,ngrid
174             zteff(ig,l)=pt(ig,l)
175        END DO
176      END DO
177      DO l=1,nlay
178        DO ig=1,ngrid
179          DO iq=1,nq
180             pqeff(ig,l,iq)=pq(ig,l,iq)
181          ENDDO
182        ENDDO
183      ENDDO
[1711]184c-----Tendencies
185      DO l=1,nlay
[1880]186        DO ig=1,ngrid
187          ztclf(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
188        ENDDO
[1711]189      ENDDO
190      DO l=1,nlay
191        DO ig=1,ngrid
192          DO iq=1,nq
[1880]193             zqclf(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
[1711]194          ENDDO
195        ENDDO
196      ENDDO
197c-----Effective temperature calculation
198      IF (CLFvarying) THEN
199         spant=3.0 ! delta T for the temprature distribution
200         mincloud=0.1 ! min cloudfrac when there is ice
201         IF (flagcloud) THEN
202             write(*,*) "Delta T", spant
203             write(*,*) "mincloud", mincloud
204             flagcloud=.false.
205         END IF
[1880]206         !CALL watersat(ngrid*nlay,ztclf,pplay,zqsat) !MV17: we dont need zqsat in the CLFvarying scheme
207         zqvap=zqclf(:,:,igcm_h2o_vap)
208         zqice=zqclf(:,:,igcm_h2o_ice)
[1711]209         CALL tcondwater(ngrid*nlay,pplay,zqvap+zqice,tcond)
210         DO l=1,nlay
211           DO ig=1,ngrid
[1880]212              zdelt=spant !MAX(spant*ztclf(ig,l),1.e-12), now totally in K. Fixed width
213              IF (tcond(ig,l) .ge. (ztclf(ig,l)+zdelt)) THEN
214                 zteff(ig,l)=ztclf(ig,l)
[1711]215                 cloudfrac(ig,l)=1.
[1880]216              ELSE IF (tcond(ig,l) .le. (ztclf(ig,l)-zdelt)) THEN
217                 zteff(ig,l)=ztclf(ig,l)-zdelt
[1711]218                 cloudfrac(ig,l)=mincloud
219              ELSE
[1880]220                 cloudfrac(ig,l)=(tcond(ig,l)-ztclf(ig,l)+zdelt)/
[1711]221     &                           (2.0*zdelt)
[1880]222                 zteff(ig,l)=(tcond(ig,l)+ztclf(ig,l)-zdelt)/2.
[1711]223              END IF
224              zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep
[1880]225              IF (cloudfrac(ig,l).le.mincloud) THEN !MV17: replaced .le.0 by .le.mincloud
[1711]226                 cloudfrac(ig,l)=mincloud
227              ELSE IF (cloudfrac(ig,l).gt.1) THEN
228                 cloudfrac(ig,l)=1.
229              END IF
230           ENDDO
231         ENDDO
232c-----Calculation of the total cloud coverage of the column
233         DO ig=1,ngrid
234            totcloudfrac(ig) = 0.
235            norm=0.
236            DO l=1,nlay
237               ponder=zqice(ig,l)*(pplev(ig,l) - pplev(ig,l+1))
238               totcloudfrac(ig) = totcloudfrac(ig)
239     &                   + cloudfrac(ig,l)*ponder
240               norm=norm+ponder
241            ENDDO
242            totcloudfrac(ig)=MAX(totcloudfrac(ig)/norm,1.e-12) ! min value if NaNs
243         ENDDO
[1880]244c-----Effective tracers quantities in the cloud fraction
245         IF (microphys) THEN
246            pqeff(:,:,igcm_ccn_mass)=pq(:,:,igcm_ccn_mass)/
247     &                              cloudfrac(:,:)
248            pqeff(:,:,igcm_ccn_number)=pq(:,:,igcm_ccn_number)/
249     &                              cloudfrac(:,:)
250         END IF ! end if (microphys)
251         pqeff(:,:,igcm_h2o_ice)=pq(:,:,igcm_h2o_ice)/
252     &                           cloudfrac(:,:)
[1711]253      END IF ! end if (CLFvarying)
[1880]254      !! MV17: CLFvarying outputs
255      CALL WRITEDIAGFI(ngrid,'pqeffice','pqeffice',
256     &             'kg/kg',3,pqeff(:,:,igcm_h2o_ice))
257      CALL WRITEDIAGFI(ngrid,'zteff','zteff',
258     &             'K',3,zteff(:,:))
259      CALL WRITEDIAGFI(ngrid,'tcond','tcond',
260     &             'K',3,tcond(:,:))
261      CALL WRITEDIAGFI(ngrid,'cloudfrac','cloudfrac',
262     &             'K',3,cloudfrac(:,:))
[633]263c------------------------------------------------------------------
[1880]264c Time subsampling for microphysics
265c------------------------------------------------------------------
[1711]266      rhocloud(1:ngrid,1:nlay) = rho_dust
[633]267      DO microstep=1,imicro
[522]268     
[633]269c-------------------------------------------------------------------
270c   1.  Tendencies:
271c------------------
[38]272
[633]273
274c------ Temperature tendency subpdt
275        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
276        ! If imicro=1 subpdt is the same as pdt
277        DO l=1,nlay
278          DO ig=1,ngrid
279             subpdt(ig,l) = subpdt(ig,l)
280     &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
281          ENDDO
282        ENDDO
283c------ Tracers tendencies subpdq
284c------ At each micro timestep we add pdq in order to have a stepped entry
285        IF (microphys) THEN
286          DO l=1,nlay
287            DO ig=1,ngrid
288              subpdq(ig,l,igcm_dust_mass) =
289     &            subpdq(ig,l,igcm_dust_mass)
290     &          + pdq(ig,l,igcm_dust_mass)
291              subpdq(ig,l,igcm_dust_number) =
292     &            subpdq(ig,l,igcm_dust_number)
293     &          + pdq(ig,l,igcm_dust_number)
294              subpdq(ig,l,igcm_ccn_mass) =
295     &            subpdq(ig,l,igcm_ccn_mass)
296     &          + pdq(ig,l,igcm_ccn_mass)
297              subpdq(ig,l,igcm_ccn_number) =
298     &            subpdq(ig,l,igcm_ccn_number)
299     &          + pdq(ig,l,igcm_ccn_number)
300            ENDDO
301          ENDDO
302        ENDIF
303        DO l=1,nlay
304          DO ig=1,ngrid
305            subpdq(ig,l,igcm_h2o_ice) =
306     &          subpdq(ig,l,igcm_h2o_ice)
307     &        + pdq(ig,l,igcm_h2o_ice)
308            subpdq(ig,l,igcm_h2o_vap) =
309     &          subpdq(ig,l,igcm_h2o_vap)
310     &        + pdq(ig,l,igcm_h2o_vap)
311          ENDDO
[1880]312        ENDDO     
[633]313       
314c-------------------------------------------------------------------
315c   2.  Main call to the different cloud schemes:
316c------------------------------------------------
317        IF (microphys) THEN
318           CALL improvedclouds(ngrid,nlay,microtimestep,
[1711]319     &             pplay,zteff,subpdt,
320     &             pqeff,subpdq,subpdqcloud,subpdtcloud,
[633]321     &             nq,tauscaling)
322
323        ELSE
324           CALL simpleclouds(ngrid,nlay,microtimestep,
[1711]325     &             pplay,pzlay,zteff,subpdt,
326     &             pqeff,subpdq,subpdqcloud,subpdtcloud,
[645]327     &             nq,tau,rice)
[633]328        ENDIF
329
330c-------------------------------------------------------------------
331c   3.  Updating tendencies after cloud scheme:
332c-----------------------------------------------
333
334        IF (microphys) THEN
335          DO l=1,nlay
336            DO ig=1,ngrid
337              subpdq(ig,l,igcm_dust_mass) =
338     &            subpdq(ig,l,igcm_dust_mass)
339     &          + subpdqcloud(ig,l,igcm_dust_mass)
340              subpdq(ig,l,igcm_dust_number) =
341     &            subpdq(ig,l,igcm_dust_number)
342     &          + subpdqcloud(ig,l,igcm_dust_number)
343              subpdq(ig,l,igcm_ccn_mass) =
344     &            subpdq(ig,l,igcm_ccn_mass)
345     &          + subpdqcloud(ig,l,igcm_ccn_mass)
346              subpdq(ig,l,igcm_ccn_number) =
347     &            subpdq(ig,l,igcm_ccn_number)
348     &          + subpdqcloud(ig,l,igcm_ccn_number)
349            ENDDO
350          ENDDO
351        ENDIF
352        DO l=1,nlay
353          DO ig=1,ngrid
354            subpdq(ig,l,igcm_h2o_ice) =
355     &          subpdq(ig,l,igcm_h2o_ice)
356     &        + subpdqcloud(ig,l,igcm_h2o_ice)
357            subpdq(ig,l,igcm_h2o_vap) =
358     &          subpdq(ig,l,igcm_h2o_vap)
359     &        + subpdqcloud(ig,l,igcm_h2o_vap)
360          ENDDO
361        ENDDO
[882]362       
363        IF (activice) THEN
364          DO l=1,nlay
365            DO ig=1,ngrid
366              subpdt(ig,l) =
367     &            subpdt(ig,l) + subpdtcloud(ig,l)
368            ENDDO
369          ENDDO
370        ENDIF
[633]371     
372 
373      ENDDO ! of DO microstep=1,imicro
374     
375c-------------------------------------------------------------------
376c   6.  Compute final tendencies after time loop:
377c------------------------------------------------
378
379c------ Temperature tendency pdtcloud
380       DO l=1,nlay
381         DO ig=1,ngrid
382             pdtcloud(ig,l) =
383     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
384          ENDDO
385       ENDDO
[740]386       
[633]387c------ Tracers tendencies pdqcloud
[703]388       DO l=1,nlay
[633]389         DO ig=1,ngrid
[703]390            pdqcloud(ig,l,igcm_h2o_ice) =
391     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro)
392     &       - pdq(ig,l,igcm_h2o_ice)
393            pdqcloud(ig,l,igcm_h2o_vap) =
394     &        subpdq(ig,l,igcm_h2o_vap)/real(imicro)
395     &       - pdq(ig,l,igcm_h2o_vap)
[740]396         ENDDO
397       ENDDO
398       
399       IF(microphys) THEN
400        DO l=1,nlay
401         DO ig=1,ngrid
[703]402            pdqcloud(ig,l,igcm_ccn_mass) =
403     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro)
404     &       - pdq(ig,l,igcm_ccn_mass)
405            pdqcloud(ig,l,igcm_ccn_number) =
406     &        subpdq(ig,l,igcm_ccn_number)/real(imicro)
407     &       - pdq(ig,l,igcm_ccn_number)
[633]408         ENDDO
[740]409        ENDDO
410       ENDIF
411       
412       IF(scavenging) THEN
413        DO l=1,nlay
414         DO ig=1,ngrid
415            pdqcloud(ig,l,igcm_dust_mass) =
416     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
417     &       - pdq(ig,l,igcm_dust_mass)
418            pdqcloud(ig,l,igcm_dust_number) =
419     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
420     &       - pdq(ig,l,igcm_dust_number)
421         ENDDO
422        ENDDO
423       ENDIF
[633]424
425c------- Due to stepped entry, other processes tendencies can add up to negative values
426c------- Therefore, enforce positive values and conserve mass
427       IF(microphys) THEN
428        DO l=1,nlay
429         DO ig=1,ngrid
[654]430          IF ((pq(ig,l,igcm_ccn_number) +
[633]431     &      ptimestep* (pdq(ig,l,igcm_ccn_number) +
[654]432     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
433     &   .or. (pq(ig,l,igcm_ccn_mass) +
434     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) +
435     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
[633]436         pdqcloud(ig,l,igcm_ccn_number) =
437     &     - pq(ig,l,igcm_ccn_number)/ptimestep
[654]438     &     - pdq(ig,l,igcm_ccn_number) + 1.
[633]439         pdqcloud(ig,l,igcm_dust_number) = 
440     &     -pdqcloud(ig,l,igcm_ccn_number)
441         pdqcloud(ig,l,igcm_ccn_mass) =
442     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
[654]443     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
[633]444         pdqcloud(ig,l,igcm_dust_mass) =
445     &     -pdqcloud(ig,l,igcm_ccn_mass)
446          ENDIF
447         ENDDO
448        ENDDO
449       ENDIF
450
[740]451       IF(scavenging) THEN
[633]452        DO l=1,nlay
453         DO ig=1,ngrid
[740]454          IF ((pq(ig,l,igcm_dust_number) +
455     &      ptimestep* (pdq(ig,l,igcm_dust_number) +
456     &        pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
457     &   .or. (pq(ig,l,igcm_dust_mass) +
458     &      ptimestep* (pdq(ig,l,igcm_dust_mass) +
459     &        pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
460         pdqcloud(ig,l,igcm_dust_number) =
461     &     - pq(ig,l,igcm_dust_number)/ptimestep
462     &     - pdq(ig,l,igcm_dust_number) + 1.
463         pdqcloud(ig,l,igcm_ccn_number) = 
464     &     -pdqcloud(ig,l,igcm_dust_number)
465         pdqcloud(ig,l,igcm_dust_mass) =
466     &     - pq(ig,l,igcm_dust_mass)/ptimestep
467     &     - pdq(ig,l,igcm_dust_mass) + 1.e-20
468         pdqcloud(ig,l,igcm_ccn_mass) =
469     &     -pdqcloud(ig,l,igcm_dust_mass)
470          ENDIF
471         ENDDO
472        ENDDO
473       ENDIF
474
475        DO l=1,nlay
476         DO ig=1,ngrid
[633]477          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
478     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
479     &       .le. 1.e-8) THEN
480           pdqcloud(ig,l,igcm_h2o_ice) =
481     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
482           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
483          ENDIF
484          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
485     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
486     &       .le. 1.e-8) THEN
487           pdqcloud(ig,l,igcm_h2o_vap) =
488     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
489           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
490          ENDIF
491         ENDDO
492        ENDDO
493
494
495c------Update the ice and dust particle size "rice" for output or photochemistry
496c------Only rsedcloud is used for the water cycle
497
498      IF(scavenging) THEN
499        DO l=1, nlay
500         DO ig=1,ngrid
501
[740]502        call updaterdust(
503     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
504     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
505     &    pdqcloud(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
506     &    pq(ig,l,igcm_dust_number) +                 ! dust number
507     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
508     &    pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
509     &    rdust(ig,l))
[633]510
511         ENDDO
512        ENDDO
[740]513      ENDIF
[1467]514
[740]515      IF(microphys) THEN
[1467]516
517       ! In case one does not want to allow supersatured water when using microphysics.
518       ! Not done by default.
519       IF(.not.supersat) THEN     
520        zt  = pt + (pdt+pdtcloud)*ptimestep
521        call watersat(ngrid*nlay,zt,pplay,zqsat)
522        DO l=1, nlay
523         DO ig=1,ngrid
524          IF (pq(ig,l,igcm_h2o_vap)
525     &      + (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
526     &      * ptimestep .ge. zqsat(ig,l)) THEN
527             pdqcloud(ig,l,igcm_h2o_vap) =
528     &         (zqsat(ig,l) - pq(ig,l,igcm_h2o_vap))/ptimestep
529     &        - pdq(ig,l,igcm_h2o_vap)
530             pdqcloud(ig,l,igcm_h2o_ice) =
531     &         -pdqcloud(ig,l,igcm_h2o_vap)
532             ! no need to correct ccn_number, updaterad can handle this properly.
533          ENDIF
534         ENDDO
535        ENDDO       
536       ENDIF
[740]537       
538       DO l=1, nlay
539         DO ig=1,ngrid
540
541        call updaterice_micro(
542     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
543     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
544     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
545     &    pq(ig,l,igcm_ccn_mass) +                   ! ccn mass
546     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
547     &    pdqcloud(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
548     &    pq(ig,l,igcm_ccn_number) +                 ! ccn number
549     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
550     &    pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
551     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
552         
[645]553         ENDDO
[740]554       ENDDO
[645]555       
[740]556      ELSE ! no microphys
557       
[645]558        DO l=1,nlay
559          DO ig=1,ngrid
[740]560         
561        call updaterice_typ(
562     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
563     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
564     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
[746]565     &    tau(ig,1),pzlay(ig,l),rice(ig,l))
[740]566
[633]567          ENDDO
[740]568         ENDDO
[633]569       
[740]570       ENDIF ! of IF(microphys)
[633]571     
[740]572     
573     
[358]574c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
575c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
576c     Then that should not affect the ice particle radius
[1047]577      do ig=1,ngrid
[358]578        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
579          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
580     &    rice(ig,2)=rice(ig,3)
581          rice(ig,1)=rice(ig,2)
582        end if
583      end do
[740]584       
585       
586       DO l=1,nlay
587         DO ig=1,ngrid
588           rsedcloud(ig,l)=max(rice(ig,l)*
589     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
590     &                    rdust(ig,l))
591!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
592         ENDDO
593       ENDDO
594       
595! used for rad. transfer calculations
596! nuice is constant because a lognormal distribution is prescribed
597      nuice(1:ngrid,1:nlay)=nuice_ref
[38]598
[1711]599c------Update tendencies for sub-grid water ice clouds
600      IF (CLFvarying) THEN
601        DO ig=1,ngrid
602          DO l=1,nlay
603            pdqcloud(ig,l,igcm_dust_mass)=pdqcloud(ig,l,igcm_dust_mass)
604     &          *cloudfrac(ig,l)
605            pdqcloud(ig,l,igcm_ccn_mass)=pdqcloud(ig,l,igcm_ccn_mass)
606     &          *cloudfrac(ig,l)
607            pdqcloud(ig,l,igcm_dust_number)=pdqcloud(ig,l,
608     &           igcm_dust_number) *cloudfrac(ig,l)
609            pdqcloud(ig,l,igcm_ccn_number)=pdqcloud(ig,l,
610     &           igcm_ccn_number) *cloudfrac(ig,l)
611            pdqcloud(ig,l,igcm_h2o_vap)=pdqcloud(ig,l,
612     &           igcm_h2o_vap) *cloudfrac(ig,l)
613            pdqcloud(ig,l,igcm_h2o_ice)=pdqcloud(ig,l,
614     &           igcm_h2o_ice) *cloudfrac(ig,l)
615          ENDDO
616        ENDDO   
617        pdtcloud(:,:)=pdtcloud(:,:)*cloudfrac(:,:)
618      ENDIF
[1758]619#ifndef MESOSCALE
[1711]620c=======================================================================
621      call WRITEDIAGFI(ngrid,"pdqice2","pdqcloudice apres microphysique"
622     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,igcm_h2o_ice))
623      call WRITEDIAGFI(ngrid,"pdqvap2","pdqcloudvap apres microphysique"
624     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
625     &      igcm_h2o_vap))
626      call WRITEDIAGFI(ngrid,"pdqccn2","pdqcloudccn apres microphysique"
627     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
628     &      igcm_ccn_mass))
629      call WRITEDIAGFI(ngrid,"pdqccnN2","pdqcloudccnN apres
630     &      microphysique","nb/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
631     &      igcm_ccn_number))
632      call WRITEDIAGFI(ngrid,"pdqdust2", "pdqclouddust apres
633     &      microphysique","kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
634     &      igcm_dust_mass))
635      call WRITEDIAGFI(ngrid,"pdqdustN2", "pdqclouddustN apres
636     &      microphysique","nb/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
637     &      igcm_dust_number))
[633]638c=======================================================================
[1758]639#endif
[633]640
[1711]641      END SUBROUTINE watercloud
642     
643      END MODULE watercloud_mod
Note: See TracBrowser for help on using the repository browser.