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

Last change on this file since 2609 was 2578, checked in by romain.vande, 4 years ago

First stage of implementing Open_MP in the physic.
So far it can initialyse physic and run with all routines at .FALSE.

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