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

Last change on this file since 2266 was 2162, checked in by emillour, 5 years ago

Mars GCM:
Some code cleanup (and preparing next steps):

  • Turn calchim into a module and make tendencies module variables in calchim_mod and watercloud_mod
  • Externalize in "physiq" the computation of solar zenithal angle (it should be computed at every physics timestep, regardless of iradia)

AB+EM

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