source: trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F90 @ 4063

Last change on this file since 4063 was 4044, checked in by jmauxion, 10 days ago

Mars PCM:
improvedclouds: from now on, the tendencies from previous physics are added
all-at-once rather than distributed over the microphysics substeps.
JM
File size: 20.7 KB
Line 
1MODULE improvedclouds_mod
2
3IMPLICIT NONE
4
5CONTAINS
6
7!=======================================================================
8
9SUBROUTINE improvedclouds(ngrid,nlay,ptimestep,pplay,pt,pdt,pq,pdq,nq,tauscaling,imicro,zt,zq)
10
11use updaterad, only: updaterice_micro, updaterccn
12use watersat_mod, only: watersat
13use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_mass, &
14                      igcm_dust_number, igcm_ccn_mass, igcm_ccn_number, igcm_hdo_vap, &
15                      igcm_hdo_ice,qperemin
16use conc_mod, only: mmean
17use comcstfi_h, only: pi, cpp
18use microphys_h, only: nbin_cld, rad_cld, mteta, kbz, nav, rgp
19use microphys_h, only: mco2, vo1, mh2o, mhdo, molco2, molhdo, To
20use nuclea_mod, only: nuclea
21use sig_h2o_mod, only: sig_h2o
22use growthrate_mod, only: growthrate
23use write_output_mod, only: write_output
24use callkeys_mod, only: activice, scavenging, cloud_adapt_ts, hdo, hdofrac
25
26implicit none
27
28!------------------------------------------------------------------
29!  This routine is used to form clouds when a parcel of the GCM is
30!    saturated. It includes the ability to have supersaturation, a
31!    computation of the nucleation rates, growthrates and the
32!    scavenging of dust particles by clouds.
33!  It is worth noting that the amount of dust is computed using the
34!    dust optical depth computed in aeropacity.F. That's why
35!    the variable called "tauscaling" is used to convert
36!    pq(dust_mass) and pq(dust_number), which are relative
37!    quantities, to absolute and realistic quantities stored in zq.
38!    This has to be done to convert the inputs into absolute
39!    values, but also to convert the outputs back into relative
40!    values which are then used by the sedimentation and advection
41!    schemes.
42
43!  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
44!           (October 2011)
45!           T. Navarro, debug,correction, new scheme (October-April 2011)
46!           A. Spiga, optimization (February 2012)
47!           J. Naar, adaptative subtimestep now done here (June 2023)
48!------------------------------------------------------------------
49
50!------------------------------------------------------------------
51!     Inputs/outputs:
52INTEGER, INTENT(IN) :: ngrid,nlay
53INTEGER, INTENT(IN) :: nq ! nombre de traceurs
54REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s)
55REAL, dimension(ngrid,nlay), INTENT(IN) :: pplay ! pression au milieu des couches (Pa)           
56REAL, dimension(ngrid,nlay), INTENT(IN) :: pt ! temperature at the middle of the layers (K)
57REAL, dimension(ngrid,nlay), INTENT(IN) :: pdt ! temperature tendency (K/s)
58REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pq ! tracer (kg/kg)
59REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pdq ! tracer tendency (kg/kg/s)
60REAL, dimension(ngrid), INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust
61INTEGER, INTENT(IN) :: imicro ! nb. microphy calls(retrocompatibility)
62
63REAL, dimension(ngrid,nlay,nq), INTENT(OUT) :: zq ! tracers post microphy (kg/kg)
64REAL, dimension(ngrid,nlay), INTENT(OUT) :: zt ! temperature post microphy (K)
65
66!------------------------------------------------------------------
67!     Local variables:
68LOGICAL, SAVE ::  firstcall = .true.
69!$OMP THREADPRIVATE(firstcall)
70REAL*8 :: derf ! Error function
71!external derf
72INTEGER :: ig,l,i
73REAL, dimension(ngrid,nlay) :: zqsat ! saturation
74REAL :: lw ! Latent heat of sublimation (J.kg-1)
75REAL :: cste
76REAL :: dMice ! mass of condensed ice
77REAL :: dMicetot ! JN
78REAL :: dMice_hdo ! mass of condensed HDO ice
79REAL, dimension(ngrid,nlay) :: alpha ! HDO equilibrium fractionation coefficient (Saturation=1)
80REAL, dimension(ngrid,nlay) :: alpha_c ! HDO real fractionation coefficient
81REAL*8 :: ph2o ! Water vapor partial pressure (Pa)
82REAL*8 :: satu ! Water vapor saturation ratio over ice
83REAL*8 :: Mo,No, Rn, Rm, dev2, n_derf, m_derf
84REAL*8, dimension(nbin_cld) :: n_aer ! number conc. of particle/each size bin
85REAL*8, dimension(nbin_cld) :: m_aer ! mass mixing ratio of particle/each size bin
86REAL :: dN, dM, seq
87REAL, dimension(nbin_cld) :: rate ! nucleation rate
88REAL, dimension(ngrid,nlay) :: rice ! Ice mass mean radius (m) (r_c in montmessin_2004)
89REAL, dimension(ngrid,nlay) :: rhocloud ! Cloud density (kg.m-3)
90REAL, dimension(ngrid,nlay) :: rdust ! Dust geometric mean radius (m)
91REAL :: res ! Resistance growth
92REAL :: Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
93
94! Parameters of the size discretization used by the microphysical scheme
95DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
96DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
97DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m)
98DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
99DOUBLE PRECISION :: vrat_cld ! Volume ratio
100DOUBLE PRECISION, dimension(nbin_cld+1), SAVE :: rb_cld ! boundary values of each rad_cld bin (m)
101!$OMP THREADPRIVATE(rb_cld)
102DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld ! width of each rad_cld bin (m)
103DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld ! particle volume for each bin (m3)
104REAL, SAVE :: sigma_ice ! Variance of the ice and CCN distributions
105!$OMP THREADPRIVATE(sigma_ice)
106
107!----------------------------------     
108! JN : used in subtimestep
109REAL :: microtimestep! subdivision of ptimestep (s)
110REAL, dimension(ngrid,nlay) :: subpdtcloud ! Temperature variation due to latent heat
111REAL, dimension(ngrid,nlay,nq) :: zq0 ! local initial value of tracers
112INTEGER, dimension(ngrid,nlay) :: zimicro ! Subdivision of ptimestep
113INTEGER, dimension(ngrid,nlay) :: count_micro ! Number of microphys calls
114REAL, dimension(ngrid,nlay) :: zpotcond ! maximal condensable water (previous two)
115REAL, dimension(1) :: zqsat_tmp ! maximal condensable water (previous two)
116REAL :: spenttime ! time spent in while loop
117REAL :: zdq ! used to compute adaptative timestep
118LOGICAL :: ending_ts ! Condition to end while loop
119
120!----------------------------------     
121! TESTS
122INTEGER :: countcells
123LOGICAL, SAVE :: test_flag    ! flag for test/debuging outputs
124!$OMP THREADPRIVATE(test_flag)
125REAL, dimension(ngrid,nlay) :: satubf,satuaf, res_out
126
127!------------------------------------------------------------------
128
129! AS: firstcall OK absolute
130IF (firstcall) THEN
131!=============================================================
132! 0. Definition of the size grid
133!=============================================================
134! rad_cld is the primary radius grid used for microphysics computation.
135! The grid spacing is computed assuming a constant volume ratio
136! between two consecutive bins; i.e. vrat_cld.
137! vrat_cld is determined from the boundary values of the size grid:
138! rmin_cld and rmax_cld.
139! The rb_cld array contains the boundary values of each rad_cld bin.
140! dr_cld is the width of each rad_cld bin.
141
142  ! Volume ratio between two adjacent bins
143  ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
144  ! vrat_cld = exp(vrat_cld)
145  vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
146  vrat_cld = exp(vrat_cld)
147  write(*,*) "vrat_cld", vrat_cld
148
149  rb_cld(1)  = rbmin_cld
150  rad_cld(1) = rmin_cld
151  vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
152
153  do i=1,nbin_cld-1
154    rad_cld(i+1)  = rad_cld(i) * vrat_cld**(1./3.)
155    vol_cld(i+1)  = vol_cld(i) * vrat_cld
156  enddo
157     
158  do i=1,nbin_cld
159    rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cld(i)
160    dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
161  enddo
162  rb_cld(nbin_cld+1) = rbmax_cld
163  dr_cld(nbin_cld)   = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
164
165  print*, ' '
166  print*,'Microphysics: size bin information:'
167  print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
168  print*,'-----------------------------------'
169  do i=1,nbin_cld
170    write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i), dr_cld(i)
171  enddo
172  write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
173  print*,'-----------------------------------'
174
175  ! we save that so that it is not computed at each timestep and gridpoint
176  rb_cld(:) = log(rb_cld(:))
177
178  ! Contact parameter of water ice on dust ( m=cos(theta) )
179  !  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
180  !  mteta is initialized in conf_phys
181  write(*,*) 'water_param contact parameter:', mteta
182
183  ! Volume of a water molecule (m3)
184  vo1 = mh2o / dble(rho_ice)
185  ! Variance of the ice and CCN distributions
186  sigma_ice = sqrt(log(1.+nuice_sed))
187     
188  write(*,*) 'Variance of ice & CCN distribs :', sigma_ice
189  write(*,*) 'nuice for sedimentation:', nuice_sed
190  write(*,*) 'Volume of a water molecule:', vo1
191
192  test_flag = .false.
193  firstcall=.false.
194ENDIF
195
196!=============================================================
197! 1. Initialisation
198!=============================================================
199
200res_out(:,:) = 0
201rice(:,:) = 1.e-8
202
203! Initialize the temperature and tracers with tendencies
204
205! If scavenging, add tendency for dust all-at-once
206IF (scavenging) THEN
207   zq(:,:,igcm_dust_mass) =zq(:,:,igcm_dust_mass)+pdq(:,:,igcm_dust_mass)*ptimestep
208   zq(:,:,igcm_dust_number) =zq(:,:,igcm_dust_number)+pdq(:,:,igcm_dust_number)*ptimestep
209ENDIF ! scavenging
210
211! Add tendency for ccn all-at-once
212zq(:,:,igcm_ccn_mass) = zq(:,:,igcm_ccn_mass) + pdq(:,:,igcm_ccn_mass)*ptimestep     
213zq(:,:,igcm_ccn_number) = zq(:,:,igcm_ccn_number) + pdq(:,:,igcm_ccn_number)*ptimestep
214
215! Add tendency for water all-at-once
216zq(:,:,igcm_h2o_ice) = zq(:,:,igcm_h2o_ice)+pdq(:,:,igcm_h2o_ice)*ptimestep
217zq(:,:,igcm_h2o_vap) = zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep
218
219! Add tendency for HDO (if computed) all-at-once
220IF (hdo) THEN
221   zq(:,:,igcm_hdo_ice) = zq(:,:,igcm_hdo_ice)+pdq(:,:,igcm_hdo_ice)*ptimestep
222   zq(:,:,igcm_hdo_vap) = zq(:,:,igcm_hdo_vap)+pdq(:,:,igcm_hdo_vap)*ptimestep
223ENDIF
224
225! Add tendency for temp all-at-one
226zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
227
228! Local temp tendency from clouds (due to latent heath release)
229subpdtcloud(:,:)=0
230
231! Handle very small values to prevent precision error
232WHERE( zq(:,:,:) < 1.e-30 ) zq(:,:,:) = 1.e-30
233
234!=============================================================
235! 2. Compute saturation
236!=============================================================
237
238dev2 = 1. / ( sqrt(2.) * sigma_ice )
239
240! Compute the condensable amount of water vapor or the sublimable water ice (if negative)
241call watersat(ngrid*nlay,max(1.,zt),pplay,zqsat) ! Make sure "temp+tendency" at least 1
242zpotcond=zq(:,:,igcm_h2o_vap) - zqsat
243     
244countcells = 0
245zimicro(:,:)=imicro
246count_micro(:,:)=0
247
248! Main loop over the GCM's grid
249DO l=1,nlay
250  DO ig=1,ngrid
251    ! Subtimestep : here we go
252    IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l), zimicro(ig,l))
253    spenttime = 0.
254    dMicetot= 0.
255    ending_ts=.false.
256    DO while (.not.ending_ts)
257      call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp)
258      zqsat(ig,l)=zqsat_tmp(1)
259      ! Get the partial pressure of water vapor and its saturation ratio
260      ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
261      satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
262      microtimestep=ptimestep/real(zimicro(ig,l))
263     
264      ! Initialize tracers for scavenging + hdo computations (JN)
265      zq0(ig,l,:)=zq(ig,l,:)
266
267      ! Check if we are integrating over ptimestep
268      if (spenttime+microtimestep.ge.ptimestep) then
269        microtimestep=ptimestep-spenttime
270        ! If so : last call !
271        ending_ts=.true.
272      endif! (spenttime+microtimestep.ge.ptimestep) then
273
274!=============================================================
275! 3. Nucleation
276!=============================================================
277
278      IF ( satu .ge. 1. ) THEN ! if there is condensation
279        call updaterccn(zq(ig,l,igcm_dust_mass),zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
280
281        ! Expand the dust moments into a binned distribution
282        Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   + 1.e-30
283        No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
284        Rn = rdust(ig,l)
285        Rn = -log(Rn)
286        Rm = Rn - 3. * sigma_ice*sigma_ice 
287        n_derf = derf( (rb_cld(1)+Rn) *dev2)
288        m_derf = derf( (rb_cld(1)+Rm) *dev2)
289        do i = 1, nbin_cld
290          n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
291          m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
292          n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
293          m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
294          n_aer(i) = n_aer(i) + 0.5 * No * n_derf
295          m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
296        enddo
297     
298        ! Get the rates of nucleation
299        call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
300     
301        dN = 0.
302        dM = 0.
303        do i = 1, nbin_cld
304          dN = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.)
305          dM = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.)
306        enddo
307
308        ! Update Dust particles
309        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
310        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
311        ! Update CCNs
312        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
313        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
314     
315      ENDIF ! of is satu >1
316
317!=============================================================
318! 4. Ice growth: scheme for radius evolution
319!=============================================================
320! We trigger crystal growth if and only if there is at least one nuclei (N>1).
321! Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
322! to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
323      IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
324        call updaterice_micro(zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),tauscaling(ig),rice(ig,l),rhocloud(ig,l))
325
326        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
327
328        ! saturation at equilibrium
329        ! rice should not be too small, otherwise seq value is not valid
330        seq  = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*max(rice(ig,l),1.e-7)))
331     
332        ! get resistance growth
333        call growthrate(zt(ig,l),pplay(ig,l),real(ph2o/satu),rice(ig,l),res,Dv)
334
335        res_out(ig,l) = res
336
337        ! implicit scheme of mass growth
338        ! cste here must be computed at each step
339        cste = 4*pi*rho_ice*microtimestep
340
341        dMice = (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))/(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
342
343        ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice.
344        dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
345        dMice = min(dMice,zq(ig,l,igcm_h2o_vap))
346
347        zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
348        zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
349
350        countcells = countcells + 1
351     
352        ! latent heat release
353        lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
354        subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
355     
356        ! DIFF: trend is enforce in a range, stabilize the scheme ?
357        if (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
358          subpdtcloud(ig,l)=5./microtimestep
359        endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
360        if (subpdtcloud(ig,l)*microtimestep.lt.-5.0) then
361            subpdtcloud(ig,l)=-5./microtimestep
362        endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
363
364        ! Special case of the isotope of water HDO   
365        if (hdo) then
366          ! condensation
367          if (dMice.gt.0.0) then
368            ! do we use fractionation?
369            if (hdofrac) then
370              ! Calculation of the HDO vapor coefficient
371              Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) ) * kbz * zt(ig,l) / ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo) * sqrt(1.+mhdo/mco2) )
372              ! Calculation of the fractionnation coefficient at equilibrium
373              alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
374              ! Calculation of the 'real' fractionnation coefficient
375              alpha_c(ig,l) = (alpha(ig,l)*satu)/( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
376            else
377              alpha_c(ig,l) = 1.d0
378            endif
379            if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
380              dMice_hdo = dMice*alpha_c(ig,l)*( zq0(ig,l,igcm_hdo_vap)/zq0(ig,l,igcm_h2o_vap) )
381            else
382              dMice_hdo=0.
383            endif
384          !! sublimation
385          else
386            if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
387              dMice_hdo=dMice*( zq0(ig,l,igcm_hdo_ice)/zq0(ig,l,igcm_h2o_ice) )
388            else
389              dMice_hdo=0.
390            endif
391          endif !if (dMice.gt.0.0)
392
393          dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
394          dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
395
396          zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
397          zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
398
399        endif ! if (hdo)       
400
401!=============================================================
402! 5. Dust cores released, tendancies, latent heat, etc ...
403!=============================================================
404
405        ! If all the ice particles sublimate, all the condensation
406        ! nuclei are released:
407        if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
408          ! Water
409          zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice)
410          zq(ig,l,igcm_h2o_ice) = 0.
411          if (hdo) then
412            zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + zq(ig,l,igcm_hdo_ice)
413            zq(ig,l,igcm_hdo_ice) = 0.
414          endif
415          ! Dust particles
416          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccn_mass)
417          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccn_number)
418          ! CCNs
419          zq(ig,l,igcm_ccn_mass) = 0.
420          zq(ig,l,igcm_ccn_number) = 0.
421        endif
422     
423      ELSE
424        ! Initialization of dMice when it's not computed
425        dMice=0 ! no condensation/sublimation to account for
426        subpdtcloud(ig,l)=0 ! no condensation/sublimation to account for
427      ENDIF !of if Nccn>1
428     
429      ! No more getting tendency : we increment tracers & temp directly
430
431      ! If not activice, the tendency from latent heat release is set to zero
432      IF (.not.activice) subpdtcloud(ig,l)=0.
433
434      ! Temperature change as a feedback from latent heat release
435      zt(ig,l) = zt(ig,l)+subpdtcloud(ig,l)*microtimestep
436
437      ! Prevent negative tracers ! JN
438      WHERE (zq(ig,l,:) < 1.e-30) zq(ig,l,:) = 1.e-30
439
440      IF (cloud_adapt_ts) then
441        ! Estimation of how much is actually condensing/subliming
442         dMicetot=dMicetot+abs(dMice) ! total accumulated ice formation since the beginning
443        IF (spenttime.ne.0) then
444          zdq=(dMicetot/spenttime)!*(ptimestep-spenttime)
445        ENDIF
446        zdq=abs(zdq)
447        call adapt_imicro(ptimestep,zdq,zimicro(ig,l))
448      ENDIF! (cloud_adapt_ts) then
449      ! Increment time spent in here
450      spenttime=spenttime+microtimestep
451      count_micro(ig,l)=count_micro(ig,l)+1
452    ENDDO ! while (.not. ending_ts)
453  ENDDO ! of ig loop
454ENDDO ! of nlayer loop
455
456!------ Useful outputs to check how it went
457call write_output("zpotcond","zpotcond microphysics","(kg/kg)",zpotcond(:,:))
458call write_output("count_micro","count_micro after microphysics","integer",count_micro(:,:))
459
460!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
461IF (test_flag) then
462  DO l=1,nlay
463    DO ig=1,ngrid
464      satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
465      satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
466    ENDDO
467  ENDDO
468  print*, 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation'
469ENDIF ! endif test_flag
470
471END SUBROUTINE improvedclouds
472
473SUBROUTINE adapt_imicro(ptimestep,potcond,zimicro)
474
475! Adaptative timestep for water ice clouds.
476! Works using a powerlaw to compute the minimal duration of subtimestep
477! (in s) should all the avalaible vapor (resp. ice) condenses (resp.sublimates)
478! Then, we use the instantaneous vap (ice) gradient extrapolated to the
479! rest of duration to increase subtimestep duration, for computing
480! efficiency.
481
482real,intent(in) :: ptimestep ! total duration of physics (sec)
483real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
484real :: alpha, beta ! Coefficients for power law
485real :: defstep,coef ! Default ptimestep of 7.5 mins (iphysiq=5)
486integer,intent(out) :: zimicro ! number of ptimestep division
487
488! Default ptimestep : defstep (7.5 mins)
489defstep=88775.*5./960.
490coef=ptimestep/defstep
491! Conservative coefficients :
492! alpha=1.81846504e+11
493! beta=1.54550140e+00
494alpha=1.88282793e+05 ! Latest values for high obliquity
495beta=4.57764370e-01  ! Latest values for high obliquity
496!alpha=1.72198978e+10 ! Present day Mars
497!beta=1.88473210e+00 
498zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.))
499!zimicro=2*zimicro ! Prediction times two, allow to complete year at high obliquity
500
501END SUBROUTINE adapt_imicro
502
503END MODULE improvedclouds_mod
Note: See TracBrowser for help on using the repository browser.