Changeset 626 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Apr 18, 2012, 10:17:12 AM (13 years ago)
Author:
tnavarro
Message:

New scheme for the clouds, no more sub-timestep. Clouds sedimentation is done with the dust one in callsedim.F like it was before. Added latent heat for sublimating ground ice. Bugs corrected. THIS VERSION OF THE WATER CYCLE SHOULD NOT BE USED WITH THERMALS DUE TO NEGATIVE TRACERS ISSUES.

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/callsedim.F

    r520 r626  
    206206        ENDIF !of if (scavenging)
    207207
    208 !        IF (water) THEN
    209 !         write(*,*) "correction for the shape of the ice particles ?"
    210 !         beta=0.75 ! default value
    211 !         call getin("ice_shape",beta)
    212 !         write(*,*) " ice_shape = ",beta
    213 !
    214 !          write(*,*) "water_param nueff Sedimentation:", nuice_sed
    215 !          IF (activice) THEN
    216 !            write(*,*) "water_param nueff Radiative:", nuice_ref
    217 !          ENDIF
    218 !        ENDIF
     208        IF (water) THEN
     209         write(*,*) "correction for the shape of the ice particles ?"
     210         beta=0.75 ! default value
     211         call getin("ice_shape",beta)
     212         write(*,*) " ice_shape = ",beta
     213
     214          write(*,*) "water_param nueff Sedimentation:", nuice_sed
     215          IF (activice) THEN
     216            write(*,*) "water_param nueff Radiative:", nuice_ref
     217          ENDIF
     218        ENDIF
    219219     
    220220        firstcall=.false.
    … …  
    383383          enddo ! of do ir=1,nr
    384384c -----------------------------------------------------------------
    385 c         WATER CYCLE CASE --- NOW DONE IN WATERCLOUD !!
    386 c -----------------------------------------------------------------
    387 !          else if (water.and.(iq.eq.igcm_h2o_ice)) then
    388 !            if (microphys) then
    389 !              ! water ice sedimentation
    390 !              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    391 !     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
    392 !     &        zqi(1,1,iq),wq,beta)
    393 !            else
    394 !              ! water ice sedimentation
    395 !              call newsedim(ngrid,nlay,ngrid*nlay,1,
    396 !     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
    397 !     &        zqi(1,1,iq),wq,beta)
    398 !            endif ! of if (microphys)
     385c         WATER CYCLE CASE
     386c -----------------------------------------------------------------
     387c          else if (water.and.(iq.eq.igcm_h2o_ice)) then
     388           else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
     389     &       .or. (iq .eq. igcm_h2o_ice)) then
     390            if (microphys) then
     391              ! water ice sedimentation
     392              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
     393     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
     394     &        zqi(1,1,iq),wq,beta)
     395            else
     396              ! water ice sedimentation
     397              call newsedim(ngrid,nlay,ngrid*nlay,1,
     398     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
     399     &        zqi(1,1,iq),wq,beta)
     400            endif ! of if (microphys)
    399401c           Tendencies
    400402c           ~~~~~~~~~~
    401 !            do ig=1,ngrid
    402 !              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
    403 !            end do
     403            do ig=1,ngrid
     404              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
     405            end do
    404406c -----------------------------------------------------------------
    405407c         GENERAL CASE
    406408c -----------------------------------------------------------------
    407 c          else
    408            else if ((iq .ne. iccn_mass) .and. (iq .ne. iccn_number)
    409      &       .and. (iq .ne. igcm_h2o_ice)) then ! because it is done in watercloud
     409          else
    410410            call newsedim(ngrid,nlay,1,1,ptimestep,
    411411     &      pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
    … …  
    445445c     Update the ice particle size "rice"
    446446c     -------------------------------------
    447 !      IF(scavenging) THEN
    448 !        DO l = 1, nlay
    449 !          DO ig=1,ngrid
    450 !            Mo   = zqi(ig,l,igcm_h2o_ice) +
    451 !     &             zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
    452 !            No   = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
    453 !            rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
    454 !     &                       +zqi(ig,l,iccn_mass)* tauscaling(ig)
    455 !     &                        / Mo * rho_dust
    456 !            rhocloud(ig,l) =
    457 !     &         min(max(rhocloud(ig,l),rho_ice),rho_dust)
    458 !            rice(ig,l) =
    459 !     &        ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
    460 !           if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8
    461 c           print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
    462 c           print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
    463 !           
    464 !          ENDDO
    465 !        ENDDO
    466 !      ELSE
    467 !        DO l = 1, nlay
    468 !          DO ig=1,ngrid
    469 !            ccntyp =
    470 !     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.)
    471 !            ccntyp = ccntyp /ccn_factor
    472 !            rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
    473 !     &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
    474 !     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
    475 !          ENDDO
    476 !        ENDDO
    477 !      ENDIF
     447      IF(scavenging) THEN
     448        DO l = 1, nlay
     449          DO ig=1,ngrid
     450            Mo   = zqi(ig,l,igcm_h2o_ice) +
     451     &             zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
     452            No   = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
     453            rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
     454     &                       +zqi(ig,l,iccn_mass)* tauscaling(ig)
     455     &                        / Mo * rho_dust
     456            rhocloud(ig,l) =
     457     &         min(max(rhocloud(ig,l),rho_ice),rho_dust)
     458            rice(ig,l) =
     459     &        ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
     460           if ((Mo.lt.1.e-15) .or. (No.le.1)) rice(ig,l) = 1.e-8
     461!           print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
     462!           print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
     463           
     464          ENDDO
     465        ENDDO
     466      ELSE
     467        DO l = 1, nlay
     468          DO ig=1,ngrid
     469            ccntyp =
     470     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.)
     471            ccntyp = ccntyp /ccn_factor
     472            rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
     473     &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
     474     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
     475          ENDDO
     476        ENDDO
     477      ENDIF
    478478
    479479      RETURN

  • trunk/LMDZ.MARS/libf/phymars/growthrate.F

    r530 r626  
    1       subroutine growthrate(timestep,temp,pmid,ph2o,psat,
    2      &                      seq,rcrystal,growth)
     1      subroutine growthrate(temp,pmid,psat,seq,rcrystal,coeff1,coeff2)
    32
    43      IMPLICIT NONE
    … …  
    109c     Authors: F. Montmessin
    1110c       Adapted for the LMD/GCM by J.-B. Madeleine (October 2011)
     11c       Use of resistances in the analytical function
     12c            instead of growth rate - T. Navarro (2012)
    1213c     
    1314c=======================================================================
    … …  
    2728c   ----------
    2829
    29       REAL timestep
     30c     Input
    3031      REAL temp     ! temperature in the middle of the layer (K)
    3132      REAL pmid     ! pressure in the middle of the layer (K)
    32       REAL*8 ph2o   ! water vapor partial pressure (Pa)
    3333      REAL*8 psat   ! water vapor saturation pressure (Pa)
     34      REAL*8 seq    ! Equilibrium saturation ratio
    3435      REAL rcrystal ! crystal radius before condensation (m)
    35       REAL*8 seq    ! Equilibrium saturation ratio
    36       REAL*8 growth
     36
     37c     Output
     38      REAL coeff1,coeff2     ! coefficients for the analytical relation between time and radius
     39
    3740
    3841c   local:
    … …  
    4346      REAL*8 afactor,Dv,lambda       ! Intermediate computations for growth rate
    4447      REAL*8 Rk,Rd
     48     
     49     
    4550
    4651c-----------------------------------------------------------------------
    … …  
    7378     &   ( pi * pmid * (molco2+molh2o)*(molco2+molh2o)
    7479     &        * sqrt(1.+mh2o/mco2) )
    75 
     80     
    7681      knudsen = temp / pmid * afactor / rcrystal
    7782      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen)
    78       Dv      = Dv / (1. + lambda * knudsen)
     83     
     84c      Dv is not corrected. Instead, we use below coefficients coeff1, coeff2
     85c      Dv      = Dv / (1. + lambda * knudsen)
    7986
    8087c     - Compute Rk
    … …  
    8289c     - Compute Rd
    8390      Rd = rgp * temp *rho_ice / (Dv*psat*mh2o)
    84 
     91     
     92     
     93      coeff1 = real(Rk + Rd*(1. + lambda * knudsen))
     94      coeff2 = real(Rk + Rd*(1. - lambda * knudsen))
     95     
     96c Below are growth rate used for other schemes
    8597c     - Compute growth=rdr/dt, then r(t+1)= sqrt(r(t)**2.+2.*growth*dt)
    86       growth = 1. / (seq*Rk+Rd)
    87 c       growth = (ph2o/psat-seq) / (seq*Rk+Rd)
     98c      growth = 1. / (seq*Rk+Rd)
     99c      growth = (ph2o/psat-seq) / (seq*Rk+Rd)
    88100c      rf   = sqrt( max( r**2.+2.*growth*timestep , 0. ) )
    89101c      dr   = rf-r

  • trunk/LMDZ.MARS/libf/phymars/improvedclouds.F

    r541 r626  
    11      subroutine improvedclouds(ngrid,nlay,ptimestep,
    2      &             pplev,pplay,zlev,pt,pdt,
     2     &             pplev,pplay,pt,pdt,
    33     &             pq,pdq,pdqcloud,pdtcloud,
    44     &             nq,tauscaling,rdust,rice,nuice,
    … …  
    2121c    values which are then used by the sedimentation and advection
    2222c    schemes.
     23c  A word about the radius growth ...
     24c  A word about nucleation and ice growth strategies ...
    2325
    2426c  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
    2527c           (October 2011)
    26 c           T. Navarro, debug & correction (October-December 2011)
     28c           T. Navarro, debug,correction, new scheme (October-April 2011)
    2729c           A. Spiga, optimization (February 2012)
    2830c------------------------------------------------------------------
    … …  
    4345      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    4446      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    45       REAL zlev(ngrid,nlay+1)    ! altitude at layer boundaries
    4647           
    4748      REAL pt(ngrid,nlay)        ! temperature at the middle of the
    … …  
    9596      REAL*8 Mo,No
    9697      REAL*8 dN,dM,newvap
    97       REAL*8 Rn, Rm, dev2, yeah, yeahn, yeahm
     98      REAL*8 Rn, Rm, dev2, yeah, n_derf, m_derf
    9899      REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
    99100      REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
    100101      REAL*8 rate(nbin_cld)  ! nucleation rate
    101       REAL*8 up,dwn,Ctot,gr,seq
     102      !REAL*8 up,dwn,Ctot,gr
     103      REAl*8 seq
    102104      REAL*8 sig      ! Water-ice/air surface tension  (N.m)
    103105      EXTERNAL sig
    … …  
    119121      SAVE sigma_ice
    120122     
     123      REAL tdicho, tmax, rmin, rmax, req, rdicho
     124      REAL coeff0, coeff1, coeff2
     125      REAL error_out(ngridmx,nlayermx)
     126      REAL error2d(ngridmx)
     127     
     128      REAL var1,var2,var3
     129     
     130      REAL rccn, epsilon
     131     
    121132c----------------------------------     
    122 c some outputs for 1D -- TEMPORARY
     133c some outputs for 1D -- TESTS
    123134      REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
    124135      REAL dN_out(ngridmx,nlayermx)  ! mass variation for output
    125136      REAL dM_out(ngridmx,nlayermx)  ! number variation for output
    126       REAL dM_col(ngridmx)           ! total mass condensed in column
    127       REAL dN_col(ngridmx)           ! total mass condensed in column
    128137      REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
    129       REAL gr_out(ngridmx,nlayermx)   ! for 1d output
    130       REAL newvap_out(ngridmx,nlayermx)   ! for 1d output
    131       REAL Mdust_col(ngridmx)        ! total column dust mass
    132       REAL Ndust_col(ngridmx)        ! total column dust mass
    133       REAL Mccn_col(ngridmx)         ! total column ccn mass
    134       REAL Nccn_col(ngridmx)         ! total column ccn mass
    135       REAL dMice_col(ngridmx)        ! total column ice mass change
    136       REAL drice_col(ngridmx)        ! total column ice radius change
    137       REAL icetot(ngridmx)        ! total column ice mass
     138      REAL gr_out(ngridmx,nlayermx)   ! for 1d output     
    138139      REAL rice_out(ngridmx,nlayermx) ! ice radius change
    139140      REAL rate_out(ngridmx,nlayermx) ! nucleation rate
     141      REAL satubf(ngridmx,nlayermx),satuaf(ngridmx,nlayermx)
     142      REAL ccnbf(ngridmx,nlayermx),ccnaf(ngridmx,nlayermx)
    140143      INTEGER count
    141144     
    142       LOGICAL output_sca     ! scavenging outputs flag for tests
    143      
    144       output_sca = .false.
     145      LOGICAL test_flag    ! flag for test/debuging outputs
     146     
     147      test_flag = .false.
    145148c----------------------------------     
    146149c----------------------------------     
    … …  
    149152
    150153      IF (firstcall) THEN
    151 
    152 c      Definition of the size grid
    153 c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     154!=============================================================
     155! 0. Definition of the size grid
     156!=============================================================
    154157c       rad_cld is the primary radius grid used for microphysics computation.
    155158c       The grid spacing is computed assuming a constant volume ratio
    … …  
    215218      END IF
    216219
    217 c-----------------------------------------------------------------------
    218 c     1. Initialization
    219 c-----------------------------------------------------------------------
     220!=============================================================
     221! 1. Initialisation
     222!=============================================================
     223
    220224c     Initialize the tendencies
    221225      pdqcloud(1:ngrid,1:nlay,1:nq)=0
    … …  
    260264          zq(ig,l,igcm_h2o_ice)=
    261265     &      pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
    262           zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004
     266          zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),1e-30) ! FF 12/2004
    263267          zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
    264268        enddo
    265269      enddo
    266270
    267 c------------------------------------------------------------------
    268 c     Cloud microphysical scheme
    269 c------------------------------------------------------------------
    270 
    271       Cste = ptimestep * 4. * pi * rho_ice
     271!=============================================================
     272! 2. Compute saturation
     273!=============================================================
     274
    272275      dev2 = 1. / ( sqrt(2.) * sigma_ice )
     276      error_out(:,:) = 0.
     277
    273278
    274279      call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
    … …  
    283288        ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
    284289        satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
    285        
    286         IF ((satu .ge. 1)!        ) THEN             ! if we have condensation
    287      &      .or. ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)
    288      &             .ge. 2)    ) THEN    ! or sublimation
    289 
     290               
     291        IF (( satu .ge. 1. )       .or.                             ! if there is condensation
     292     &      ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.)) THEN ! or sublimation
     293
     294
     295!=============================================================
     296! 3. Nucleation
     297!=============================================================
    290298
    291299c       Expand the dust moments into a binned distribution
    292         Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
    293         No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30
     300        Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   !+ 1.e-30
     301        No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
    294302        Rn = rdust(ig,l)
    295303        Rn = -dlog(Rn)
    296304        Rm = Rn - 3. * sigma_ice*sigma_ice 
    297         yeahn = derf( (rb_cld(1)+Rn) *dev2)
    298         yeahm = derf( (rb_cld(1)+Rm) *dev2)
     305        n_derf = derf( (rb_cld(1)+Rn) *dev2)
     306        m_derf = derf( (rb_cld(1)+Rm) *dev2)
    299307        do i = 1, nbin_cld
    300           n_aer(i) = -0.5 * No * yeahn !! this ith previously computed
    301           m_aer(i) = -0.5 * Mo * yeahm !! this ith previously computed
    302           yeahn = derf( (rb_cld(i+1)+Rn) *dev2)
    303           yeahm = derf( (rb_cld(i+1)+Rm) *dev2)
    304           n_aer(i) = n_aer(i) + 0.5 * No * yeahn
    305           m_aer(i) = m_aer(i) + 0.5 * Mo * yeahm
     308          n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
     309          m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
     310          n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
     311          m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
     312          n_aer(i) = n_aer(i) + 0.5 * No * n_derf
     313          m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
    306314        enddo
    307315!!! MORE EFFICIENT COMPUTATIONALLY THAN
    … …  
    350358        rate_out(ig,l) = 0
    351359        do i = 1, nbin_cld
    352         ! schema simple
    353360          n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
    354361          m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
    … …  
    358365        enddo
    359366
    360           Mo   = zq0(ig,l,igcm_h2o_ice) +
    361      &           zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
    362           No   = zq0(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
    363           !write(*,*) "l,cloud particles,cloud mass",l, No, Mo
    364           rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
    365      &                     +zq0(ig,l,igcm_ccn_mass)* tauscaling(ig)
    366      &                      / Mo * rho_dust
    367           rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
    368           if ((Mo.lt.1.e-20) .or. (No.le.1)) then
    369               rice(ig,l) = 1.e-8
     367c       Update CCNs, can also be done after the radius growth
     368c       Dust particles
     369        zq(ig,l,igcm_dust_mass)   =
     370     &         zq(ig,l,igcm_dust_mass)   - dM/ tauscaling(ig)
     371        zq(ig,l,igcm_dust_number) =
     372     &         zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig)
     373c       CCNs
     374        zq(ig,l,igcm_ccn_mass)   =
     375     &         zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig)
     376        zq(ig,l,igcm_ccn_number) =
     377     &         zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
     378       
     379
     380!=============================================================
     381! 4. Ice growth: scheme for radius evolution
     382!=============================================================
     383
     384        Mo   = zq(ig,l,igcm_h2o_ice) +
     385     &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig)  + 1.e-30
     386        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
     387
     388        rhocloud(ig,l) =  zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
     389     &                  + zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
     390     &                  / Mo * rho_dust
     391        rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
     392         
     393c       nuclei radius 
     394        rccn  = CBRT(zq(ig,l,igcm_ccn_mass)/
     395     &                (pi*rho_dust*zq(ig,l,igcm_ccn_number)*4./3.))
     396
     397c       ice crystal radius   
     398        rice (ig,l) =
     399     &      CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) )
     400
     401c       enforce physical value : crystal cannot be smaller than its nuclei !
     402        rice(ig,l) = max(rice(ig,l), rccn)   
     403
     404c       saturation at equilibrium
     405        seq  = exp( 2.*sig(zt(ig,l))*mh2o /
     406     &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )
     407
     408
     409c       If there is more than on nuclei, we peform ice growth 
     410        var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
     411        IF (var1 .ge. -1) THEN
     412
     413
     414      if (test_flag) then
     415       print*, ' '
     416       print*, ptimestep
     417       print*, 'satu,seq', real(satu), real(seq), ig,l
     418       print*, 'dN,dM', real(dN),real(dM)           
     419       print*,'rccn', rccn
     420       print*, 'Nccn, Mccn', zq(ig,l,igcm_ccn_number)*tauscaling(ig),
     421     &  zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
     422      endif
     423     
     424c         crystal radius to reach saturation at equilibrium (i.e. satu=seq)         
     425          req = ( zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
     426     &          + zq(ig,l,igcm_ccn_mass)*tauscaling(ig)*
     427     &           rho_ice/rho_dust - seq * zqsat(ig,l))
     428     &         / ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)*
     429     &             pi*rho_ice*4./3. )
     430          req = CBRT(req)
     431         
     432c         compute ice radius growth resistances (diffusive and latent heat resistancea)       
     433          call growthrate(zt(ig,l),pplay(ig,l),
     434     &                    ph2o/satu,seq,req,coeff1,coeff2)
     435               
     436          coeff0 = -zqsat(ig,l) / (4.*pi*req*rho_ice
     437     &               * zq(ig,l,igcm_ccn_number)*tauscaling(ig))
     438               
     439c         compute tmax, the time needed to reach req           
     440          call phi(rice(ig,l),req,coeff1,coeff2,tmax)
     441         
     442      if (test_flag) then
     443          print*, 'coeffs', coeff0,coeff1,coeff2
     444          print*, 'req,tmax', req,tmax*coeff0
     445          print*, 'i,rmin,rdicho,rmax,tdicho'
     446      endif
     447
     448c         rmin is rice if r increases (satu >1) or req if it decreases (satu<1)
     449c         if req is lower than rccn (ie not enough ice to reach saturation), rmin is forced to rccn
     450          if (satu .ge. seq) then
     451            ! crystal size is increasing
     452            rmin = max(min(rice(ig,l),req),rccn)
     453            rmax = max(rice(ig,l),req)
    370454          else
    371               rice(ig,l) =
    372      &          CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
     455            !  crystal size is decreasing
     456            rmax = max(min(rice(ig,l),req),rccn)
     457            rmin = max(rice(ig,l),req)
    373458          endif
    374 c          nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
    375           seq  = exp( 2.*sig(zt(ig,l))*mh2o /
    376      &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )
    377 
    378           call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
    379      &                    ph2o,ph2o/satu,seq,rice(ig,l),gr)
    380 
    381           up  = Cste * gr * rice(ig,l) * No * seq +
    382      &            zq(ig,l,igcm_h2o_vap)
    383           dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
    384  
    385           Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
    386          
    387           newvap = min(up/dwn,Ctot)
    388  
    389           gr = gr * ( newvap/zqsat(ig,l) - seq )
    390 
    391          
    392           dMice = min( max(Cste * No * rice(ig,l) * gr,
    393      &                     -zq(ig,l,igcm_h2o_ice) ),
    394      &                 zq(ig,l,igcm_h2o_vap) )
    395      
    396 c----------- TESTS 1D output ---------
    397              satu_out(ig,l) = satu
    398              Mcon_out(ig,l) = dMice
    399              newvap_out(ig,l) = newvap
    400              gr_out(ig,l) = gr
    401              dN_out(ig,l) = dN
    402              dM_out(ig,l) = dM
    403 c-------------------------------------
    404          
    405 c       Water ice
    406           zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) +
    407      &        dMice
    408      
    409 c       Water vapor
    410           zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) -
    411      &        dMice
    412      
     459                         !rmax = min(rmax,1.e-3) ! au max on a des rayons de 1 mm  pour la dicho ...
     460          rdicho = 0.5*(rmin+rmax)
     461         
     462          ! for output
     463          var1 = rice(ig,l)
     464          var2 = rmin
     465          var3 = rmax
     466         
     467c         Given the phi function is monotonous, we perform a simple dichotomy to find the raidus at t+1       
     468          do i = 1,10 ! dichotomy loop
     469             
     470c            compute tdicho, the time needed to reach rdicho
     471             call phi(rdicho,req,coeff1,coeff2,tdicho)
     472             !print*, tdicho,tmax
     473             tdicho = coeff0*(tdicho - tmax)
     474             
     475             if (test_flag) print*, i,rmin,rdicho,rmax,tdicho
     476                       
     477             if (tdicho .ge. ptimestep) then
     478                rmax = rdicho
     479             else
     480                rmin = rdicho
     481             endif
     482             
     483             rdicho = 0.5*(rmin+rmax)
     484         
     485          enddo  ! of dichotomy loop
     486         
     487       if (test_flag) then
     488        if (abs(rdicho - rccn) .ge. 1e-15) then ! to avoid infinite values
     489          epsilon = (rmax - rmin)/(2**10)
     490          error_out(ig,l) =
     491     &    100*(epsilon**3 +3*epsilon**2*rdicho +3*epsilon*rdicho**2)
     492     &    / (rdicho**3-rccn**3)     
     493        endif
     494        print*, 'error masse glace %', error_out(ig,l)
     495        print*, 'rice,ice,vap bf',
     496     &   rice(ig,l),zq0(ig,l,igcm_h2o_ice),zq0(ig,l,igcm_h2o_vap)
     497       endif
     498         
     499c         If the initial condition is subsaturated and there is not enough ice available for sublimation
     500c         to reach equilibrium, req is neagtive. Therefore, enforce physical value.
     501          rice(ig,l) = max(rdicho,rccn)
     502         
     503!!!!! Water ice mass is computed with radius at t+1 and their current number
     504!!!!! Nccn is at t or t+1, depending on what has been done before
     505!          zq(ig,l,igcm_h2o_ice) =
     506!     &       (pi*rho_ice*zq(ig,l,igcm_ccn_number)*4./3.
     507!     &         * rdicho*rdicho*rdicho -
     508!     &         zq(ig,l,igcm_ccn_mass)*rho_ice/rho_dust)
     509!     &       * tauscaling(ig)
     510
     511!!!!! Water ice mass is computed with radius at t+1 and their number at t+1
     512!!!!! that is dirty, but this enforces the use of Nccn at t+1, whatever is done before
     513!!!!! TO BE CLEANED ONE DAY
     514          var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
     515          var2 = zq0(ig,l,igcm_ccn_mass)*tauscaling(ig)   + dM
     516          zq(ig,l,igcm_h2o_ice) =
     517     &       (pi*rho_ice*var1*4./3.
     518     &         * rdicho*rdicho*rdicho -
     519     &         var2*rho_ice/rho_dust)
     520     
     521     
     522      !!!! enforce realistic values (if the case of growth on Nccn(t) and condensation on Nccn(t+1))
     523          zq(ig,l,igcm_h2o_ice) =
     524     &     min(zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap),
     525     &     zq(ig,l,igcm_h2o_ice))
     526          zq(ig,l,igcm_h2o_ice) =
     527     &     max(1e-30,zq(ig,l,igcm_h2o_ice))
     528     
     529          zq(ig,l,igcm_h2o_vap) =
     530     &         zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
     531     &       - zq(ig,l,igcm_h2o_ice)
     532     
     533     
     534       if (test_flag) then
     535        print*, 'rice,ice,vap af',
     536     &     rice(ig,l),zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_h2o_vap)
     537        print*, 'satu bf, af',
     538     &     zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l),
     539     &     zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
     540       endif
     541     
     542     
     543!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST           
     544          if ((zq(ig,l,igcm_h2o_ice).le. -1e-8)
     545     &    .or. (zq(ig,l,igcm_h2o_vap).le. -1e-8)) then
     546           print*, 'NEGATIVE WATER'
     547           print*, 'ig,l', ig,l
     548           print*, 'satu', satu
     549           print*, 'vap, ice bf',
     550     &          zq0(ig,l,igcm_h2o_vap), zq0(ig,l,igcm_h2o_ice)
     551           print*, 'vap, ice af',
     552     &          zq(ig,l,igcm_h2o_vap), zq(ig,l,igcm_h2o_ice)
     553           print*, 'ccn N,M bf',
     554     &          zq0(ig,l,igcm_ccn_number), zq0(ig,l,igcm_ccn_mass)
     555           print*, 'ccn N,M af',
     556     &          zq(ig,l,igcm_ccn_number), zq(ig,l,igcm_ccn_mass)
     557           print*, 'tauscaling',
     558     &          tauscaling(ig)
     559           print*, 'req,rccn,rice bf,rice af',
     560     &         req,rccn,var1,rice(ig,l)
     561           print*, 'rmin, rmax', var2,var3
     562           print*, 'error_out,tdicho,ptimestep',
     563     &          error_out(ig,l),tdicho,ptimestep
     564           print*, ' '
     565          endif
     566!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST           
     567         
     568         
     569        ENDIF !of if Nccn >1
     570         
     571     
     572!=============================================================
     573! 5. Dust cores released, tendancies, latent heat, etc ...
     574!=============================================================
     575
    413576c         If all the ice particles sublimate, all the condensation
    414 c           nuclei are released:
    415           if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
    416 c           for coherence
    417             dM = 0
    418             dN = 0
    419             dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
    420             dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
    421 c           Water ice particles
     577c         nuclei are released:
     578          if (zq(ig,l,igcm_h2o_ice).le.1e-10) then
     579!           for coherence
     580!            dM = 0
     581!            dN = 0
     582!            dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
     583!            dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
     584c           Water
     585            zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
     586     &                            + zq(ig,l,igcm_h2o_ice)
    422587            zq(ig,l,igcm_h2o_ice) = 0.
    423588c           Dust particles
    424             zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) +
    425      &        zq(ig,l,igcm_ccn_mass)
    426             zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) +
    427      &        zq(ig,l,igcm_ccn_number)
     589            zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
     590     &                              + zq(ig,l,igcm_ccn_mass)
     591            zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
     592     &                                + zq(ig,l,igcm_ccn_number)
    428593c           CCNs
    429594            zq(ig,l,igcm_ccn_mass) = 0.
    430595            zq(ig,l,igcm_ccn_number) = 0.
    431596          endif
    432 
    433         dN = dN/ tauscaling(ig)
    434         dM = dM/ tauscaling(ig)
    435 c       Dust particles
    436         zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
    437         zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
    438 c       CCNs
    439         zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
    440         zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
     597         
     598
     599!        dN = dN/ tauscaling(ig)
     600!        dM = dM/ tauscaling(ig)
     601!c       Dust particles
     602!        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
     603!        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
     604!c       CCNs
     605!        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
     606!        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
    441607
    442608               
    … …  
    456622     
    457623        count = count +1
     624       
     625       
    458626        ELSE ! for coherence (rdust, rice computations etc ..)
    459627          zq(ig,l,igcm_dust_mass)   = zq0(ig,l,igcm_dust_mass)
    … …  
    463631          zq(ig,l,igcm_h2o_ice)     = zq0(ig,l,igcm_h2o_ice)
    464632          zq(ig,l,igcm_h2o_vap)     = zq0(ig,l,igcm_h2o_vap)
    465 
     633         
    466634! pour les sorties de test
    467           satu_out(ig,l) = satu
    468           Mcon_out(ig,l) = 0
    469           newvap_out(ig,l) = zq(ig,l,igcm_h2o_vap)
    470           gr_out(ig,l) = 0
    471           dN_out(ig,l) = 0
    472           dM_out(ig,l) = 0
     635!          satu_out(ig,l) = satu
     636!          gr_out(ig,l) = 0
     637!          dN_out(ig,l) = 0
     638!          dM_out(ig,l) = 0
    473639         
    474640        ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
    475641       
    476 c-----update temperature       
     642!        ccnbf(ig,l) = zq0(ig,l,igcm_ccn_number) * tauscaling(ig)
     643!        ccnaf(ig,l) = zq(ig,l,igcm_ccn_number) * tauscaling(ig)
     644!
     645!        satubf(ig,l) = zq0(ig,l,igcm_h2o_vap) / zqsat(ig,l)
     646!        satuaf(ig,l) = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
     647       
     648       
     649c-----update temperature (latent heat relase)     
    477650          lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
    478           pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
    479 c          pdtcloud(ig,l)=pdt(ig,l)-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
     651          pdtcloud(ig,l)= -pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
    480652         
    481653c----- update rice & rhocloud AFTER microphysic
    … …  
    488660          rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
    489661         
    490           rice_out(ig,l)=rice(ig,l)
    491662          if ((Mo.lt.1.e-20) .or. (No.le.1)) then
    492663              rice(ig,l) = 1.e-8
    … …  
    495666     &          CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
    496667          endif
    497           rice_out(ig,l)=rice(ig,l)-rice_out(ig,l)
    498668         
    499669          nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
    … …  
    510680          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
    511681
    512         ENDDO
    513       ENDDO
    514      
    515 !      print*, 'SATU'
    516 !      print*, satu_out(1,:)
    517      
    518 c------------------------------------------------------------------
    519 c------------------------------------------------------------------
    520 c------------------------------------------------------------------
    521 c------------------------------------------------------------------
    522 c------------------------------------------------------------------
    523 c     TESTS
    524 
    525       IF (output_sca) then
    526  
     682        ENDDO ! of ig loop
     683      ENDDO ! of nlayer loop
     684     
     685     
     686!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
     687!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
     688!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
     689      IF (test_flag) then
     690     
     691       error2d(:) = 0.
     692       DO l=1,nlay
     693       DO ig=1,ngrid
     694         error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
     695       ENDDO
     696       ENDDO
     697
    527698      print*, 'count is ',count, ' i.e. ',
    528699     &     count*100/(nlay*ngrid), '% for microphys computation'     
    529700
    530       dM_col(:)    = 0
    531       dN_col(:)    = 0
    532       Mdust_col(:) = 0
    533       Ndust_col(:) = 0
    534       Mccn_col(:)  = 0
    535       Nccn_col(:)  = 0
    536       dMice_col(:) = 0
    537       drice_col(:) = 0
    538       icetot(:)    = 0
    539       do l=1, nlay
    540         do ig=1,ngrid
    541           dM_col(ig) = dM_col(ig) +
    542      &       dM_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
    543           dN_col(ig) = dN_col(ig) +
    544      &       dN_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
    545           Mdust_col(ig) = Mdust_col(ig) +
    546      &       zq(ig,l,igcm_dust_mass)*tauscaling(ig)
    547      &       *(pplev(ig,l) - pplev(ig,l+1)) / g
    548           Ndust_col(ig) = Ndust_col(ig) +
    549      &       zq(ig,l,igcm_dust_number)*tauscaling(ig)
    550      &       *(pplev(ig,l) - pplev(ig,l+1)) / g
    551           Mccn_col(ig) = Mccn_col(ig) +
    552      &       zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
    553      &       *(pplev(ig,l) - pplev(ig,l+1)) / g
    554           Nccn_col(ig) = Nccn_col(ig) +
    555      &       zq(ig,l,igcm_ccn_number)*tauscaling(ig)
    556      &       *(pplev(ig,l) - pplev(ig,l+1)) / g
    557           dMice_col(ig) = dMice_col(ig) +
    558      &       Mcon_out(ig,l)
    559      &       *(pplev(ig,l) - pplev(ig,l+1)) / g
    560           drice_col(ig) = drice_col(ig) +
    561      &       rice_out(ig,l)*zq(ig,l,igcm_h2o_ice)
    562      &       *(pplev(ig,l) - pplev(ig,l+1)) / g
    563           icetot(ig) = icetot(ig) +
    564      &       zq(ig,l,igcm_h2o_ice)
    565      &       *(pplev(ig,l) - pplev(ig,l+1)) / g
    566         enddo ! of do ig=1,ngrid
    567       enddo ! of do l=1,nlay
    568      
    569       drice_col(ig) = drice_col(ig)/icetot(ig)
    570 
    571701      IF (ngrid.ne.1) THEN ! 3D
    572          call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
    573      &                    satu_out)
    574          call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",2,
    575      &                    dM_col)
    576          call WRITEDIAGFI(ngrid,"dN_col","dN column","N",2,
    577      &                    dN_col)
    578          call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",2,
    579      &                    Ndust_col)
    580          call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",2,
    581      &                    Mdust_col)
    582          call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",2,
    583      &                    Nccn_col)
    584          call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",2,
    585      &                    Mccn_col)
    586          call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
    587      &                    dM_out)
    588          call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
    589      &                    dN_out)
     702!         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
     703!     &                    satu_out)
     704!         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
     705!     &                    dM_out)
     706!         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
     707!     &                    dN_out)
     708         call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
     709     &                    error2d)
     710!         call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
     711!     &                    zqsat)
    590712      ENDIF
    591713
    592714      IF (ngrid.eq.1) THEN ! 1D
    593 
    594          call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
    595      &                    newvap_out)
    596          call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
    597      &                    gr_out)
    598          call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
    599      &                    rate_out)
    600          call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
    601      &                    dM_out)
    602          call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
    603      &                    dN_out)
    604          call WRITEDIAGFI(ngrid,"dicetot","ice col var","kg/m2",0,
    605      &                    dMice_col)
    606          call WRITEDIAGFI(ngrid,"dricetot","ice col var","m",0,
    607      &                    drice_col)
    608          call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
    609      &                    Mcon_out)
     715         call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
     716     &                    error_out)
     717         call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1,
     718     &                    satubf)
     719         call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1,
     720     &                    satuaf)
     721         call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1,
     722     &                    zq0(1,:,igcm_h2o_vap))
     723         call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1,
     724     &                    zq(1,:,igcm_h2o_vap))
     725         call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1,
     726     &                    zq0(1,:,igcm_h2o_ice))
     727         call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1,
     728     &                    zq(1,:,igcm_h2o_ice))
     729         call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1,
     730     &                    ccnbf)
     731         call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1,
     732     &                    ccnaf)
     733c         call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
     734c     &                    gr_out)
     735c         call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
     736c     &                    rate_out)
     737c         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
     738c     &                    dM_out)
     739c         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
     740c     &                    dN_out)
    610741         call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
    611742     &                    zqsat)
    … …  
    620751         call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
    621752     &                    rhocloud)
    622          call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",0,
    623      &                    dM_col)
    624          call WRITEDIAGFI(ngrid,"dN_col","dN column","N",0,
    625      &                    dN_col)
    626          call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",0,
    627      &                    Ndust_col)
    628          call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",0,
    629      &                    Mdust_col)
    630          call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",0,
    631      &                    Nccn_col)
    632          call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",0,
    633      &                    Mccn_col)
    634753      ENDIF
    635       ENDIF ! endif output_sca
    636 c------------------------------------------------------------------
     754     
     755      ENDIF ! endif test_flag
     756!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
     757!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
     758!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
     759
    637760      return
    638761      end
     762     
     763     
     764     
     765cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
     766cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
     767c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
     768c It is an analytical solution to the ice radius growth equation,
     769c with the approximation of a constant 'reduced' cunningham correction factor
     770c (lambda in growthrate.F) taken at radius req instead of rice   
     771cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
     772cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     773
     774      subroutine phi(rice,req,coeff1,coeff2,time)
     775     
     776      implicit none
     777     
     778      ! inputs
     779      real rice ! ice radius
     780      real req  ! ice radius at equilibirum
     781      real coeff1  ! coeff for the log
     782      real coeff2  ! coeff for the arctan
     783
     784      ! output     
     785      real time
     786     
     787      !local
     788      real var
     789     
     790      ! 1.73205 is sqrt(3)
     791     
     792      var = max(
     793     &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
     794           
     795       time =
     796     &   coeff1 *
     797     &   log( var )
     798     & + coeff2 * 1.73205 *
     799     &   atan( (2*rice+req) / (1.73205*req) )
     800     
     801      return
     802      end
     803     
     804     
     805     

  • trunk/LMDZ.MARS/libf/phymars/lwu.F

    r38 r626  
    2727c
    2828c MODIF : FF : removing the monster bug on water ice clouds 11/2010
     29c
     30c MODIF : TN : corrected bug if very big water ice clouds 04/2012
    2931c-----------------------------------------------------------------------
    3032 
    … …  
    201203c   et pourquoi pas d'abord?  hourdin@lmd.ens.fr
    202204
    203            aer_t(jl,ja,jkl)=exp(-pview*aer_a(jl,ja,jkl))
     205c TN 04/12 : if very big water ice clouds, aer_t is strictly rounded
     206c to zero in lower levels, which is a source of NaN
     207           !aer_t(jl,ja,jkl)=exp(-pview*aer_a(jl,ja,jkl))
     208           aer_t(jl,ja,jkl)=max(exp(-pview*aer_a(jl,ja,jkl)),1e-30)
     209           
    204210         
    205211          enddo
    206212        enddo
    207       enddo                 
    208 
     213      enddo     
     214     
    209215c----------------------------------------------------------------------
    210216      return

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r608 r626  
    194194      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
    195195     
    196       REAL watercapflag(ngridmx)     ! water cap flag
    197 
    198196c     Variables used by the water ice microphysical scheme:
    199197      REAL rice(ngridmx,nlayermx)    ! Water ice geometric mean radius (m)
    … …  
    986984           call watercloud(ngrid,nlayer,ptimestep,
    987985     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
    988      &                pq,pdq,zdqcloud,zdqscloud,zdtcloud,
     986     &                pq,pdq,zdqcloud,zdtcloud,
    989987     &                nq,tau,tauscaling,rdust,rice,nuice,
    990988     &                rsedcloud,rhocloud)
    … …  
    10081006     &                    pdq(ig,l,igcm_h2o_ice)+
    10091007     &                    zdqcloud(ig,l,igcm_h2o_ice)
    1010                  if (((pq(ig,l,igcm_h2o_ice) +
     1008             ! There are negative values issues with some tracers (17/04)
     1009                 if (((pq(ig,l,igcm_h2o_ice) +
    10111010     &                 pdq(ig,l,igcm_h2o_ice)*ptimestep)) .le. 0)
    1012      &           then
    1013                    pdq(ig,l,igcm_h2o_ice) =
    1014      &                 - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30
    1015                  endif
     1011     &             pdq(ig,l,igcm_h2o_ice) =               
     1012     &                 - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30 
     1013                 if (((pq(ig,l,igcm_h2o_vap) +
     1014     &                 pdq(ig,l,igcm_h2o_vap)*ptimestep)) .le. 0)
     1015     &             pdq(ig,l,igcm_h2o_vap) =               
     1016     &                 - pq(ig,l,igcm_h2o_vap)/ptimestep + 1e-30
     1017     
    10161018                 IF (scavenging) THEN
    10171019                   pdq(ig,l,igcm_ccn_mass)=
    … …  
    10521054           ENDIF ! of IF (water) THEN
    10531055
    1054 ! Increment water ice surface tracer tendency
    1055          DO ig=1,ngrid
    1056            dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
    1057      &                               zdqscloud(ig,igcm_h2o_ice)
    1058          ENDDO
    1059          
    10601056         END IF  ! of IF (water)
    10611057
    … …  
    18231819     &                       'surface h2o_ice',
    18241820     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
    1825 
    1826 c            if (caps) then
    1827 c             do ig=1,ngridmx
    1828 c                if (watercaptag(ig)) watercapflag(ig) = 1
    1829 c             enddo
    1830 c            CALL WRITEDIAGFI(ngridmx,'watercaptag',
    1831 c     &                         'Ice water caps',
    1832 c     &                         '',2,watercapflag)
    1833 cs            endif
    18341821c           CALL WRITEDIAGFI(ngridmx,'albedo',
    18351822c    &                         'albedo',

  • trunk/LMDZ.MARS/libf/phymars/updatereffrad.F

    r420 r626  
    6969      REAL, PARAMETER :: ccn0 = 1.3E8
    7070
    71       LOGICAL firstcall
    72       DATA firstcall/.true./
    73       SAVE firstcall
     71c      LOGICAL firstcall
     72c      DATA firstcall/.true./
     73c      SAVE firstcall
    7474
    7575      REAL CBRT
    … …  
    8383c==================================================================
    8484
    85       IF (firstcall) THEN
     85c      IF (firstcall) THEN
    8686c       At firstcall, rdust and rice are not known; therefore
    8787c         they need to be computed below.
     88
     89c      Correction TN 17/04: rdust and rice must be updated at all steps,
     90c      otherwise it is a possible source of bugs
    8891
    8992c       1.1 Dust particles
    … …  
    120123          ENDDO
    121124        ENDIF ! of if (water.AND.activice)
    122         firstcall = .false.
    123       ENDIF ! of if firstcall
     125       
     126c        firstcall = .false.
     127c      ENDIF ! of if firstcall
    124128
    125129c==================================================================

  • trunk/LMDZ.MARS/libf/phymars/vdifc.F

    r557 r626  
    3939#include "comgeomfi.h"
    4040#include "tracer.h"
     41#include "microphys.h"
    4142
    4243c
    … …  
    99100      LOGICAL firstcall
    100101      SAVE firstcall
     102
     103      REAL lw
    101104
    102105c     variable added for CO2 condensation:
    … …  
    796799c             Starting upward calculations for water :
    797800               zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
     801c             Take into account H2o latent heat in surface energy budget         
     802               lw = (2834.3-0.28*(ptsrf(ig)-To)
     803     &              -0.004*(ptsrf(ig)-To)**2)*1.e+3
     804               pdtsrf(ig) = pdtsrf(ig)
     805     &                    + pdqsdif(ig,igcm_h2o_ice)*lw/pcapcal(ig)
    798806            ENDDO ! of DO ig=1,ngrid
    799807          ELSE
    … …  
    806814          DO ilay=2,nlay
    807815             DO ig=1,ngrid
    808                 zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
     816               zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
    809817             ENDDO
    810818          ENDDO

  • trunk/LMDZ.MARS/libf/phymars/watercloud.F

    r576 r626  
    1        SUBROUTINE watercloud(ngrid,nlay,ptimestep,
     1       SUBROUTINE watercloud(ngrid,nlay, ptimestep,
    22     &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
    3      &                pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     3     &                pq,pdq,pdqcloud,pdtcloud,
    44     &                nq,tau,tauscaling,rdust,rice,nuice,
    55     &                rsedcloud,rhocloud)
    6 ! to use  'getin'
    7       USE ioipsl_getincom
    86      IMPLICIT NONE
    97
    … …  
    1513c    - An improved microphysical scheme (see improvedclouds.F)
    1614c
    17 c  There is a time loop specific to cloud formation and sedimentation
    18 c  due to timescales smaller than the GCM integration timestep.
    19 c
    2015c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
    2116c           J.-B. Madeleine, Thomas Navarro
    2217c
    23 c  2004 - 2012
     18c  2004 - Oct. 2011
    2419c=======================================================================
    2520
    … …  
    4035
    4136      INTEGER ngrid,nlay
    42       INTEGER nq                 ! nombre de traceurs
     37      integer nq                 ! nombre de traceurs
    4338      REAL ptimestep             ! pas de temps physique (s)
    4439      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    4540      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    46       REAL pdpsrf(ngrid)         ! tendence surf pressure
     41      REAL pdpsrf(ngrid)         ! tendance surf pressure
    4742      REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
    4843      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
    4944      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
    50       REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
     45      REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.
    5146
    5247      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
    53       real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
     48      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
    5449
    55       REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
    56       REAL tauscaling(ngridmx)   ! Convertion factor for dust amount
    57       real rdust(ngridmx,nlay)  ! Dust geometric mean radius (m)
     50      REAL tau(ngridmx,naerkind)   ! Column dust optical depth at each point
     51      REAL tauscaling(ngridmx)     ! Convertion factor for dust amount
     52      real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
    5853
    5954c   Outputs:
    6055c   -------
    6156
    62       real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
    63       real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
    64       REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
    65                                    ! a la chaleur latente
     57      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
     58      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
     59                                   !   a la chaleur latente
    6660
    6761      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
    … …  
    6963      REAL nuice(ngrid,nlay)   ! Estimated effective variance
    7064                               !   of the size distribution
    71       real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
    72       real rhocloud(ngridmx,nlay)  ! Cloud density (kg.m-3)
     65      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
     66      real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
    7367
    7468c   local:
    7569c   ------
    7670
    77       ! for sedimentation
    78       REAL zq(ngridmx,nlay,nqmx)    ! local value of tracers
    79       real masse (ngridmx,nlay)     ! Layer mass (kg.m-2)
    80       real epaisseur (ngridmx,nlay) ! Layer thickness (m)
    81       real wq(ngridmx,nlay+1)       ! displaced tracer mass (kg.m-2)
    82      
    83       ! for ice radius computation
    84       REAL Mo,No
    85       REAl ccntyp
    86       real beta      ! correction for the shape of the ice particles (cf. newsedim)
    87       save beta
    88      
    89       ! for time loop
    90       INTEGER microstep  ! time subsampling step variable
    91       INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
    92       SAVE imicro
    93       REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
    94       SAVE microtimestep
    95      
    96       ! tendency given by clouds (inside the micro loop)
    97       REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
    98       REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
    99 
    100       ! global tendency (clouds+sedim+physics)
    101       REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
    102       REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
    103      
    104       REAL CBRT
    105       EXTERNAL CBRT
    106 
    107      
    108       INTEGER iq,ig,l
     71      INTEGER ig,l
    10972      LOGICAL,SAVE :: firstcall=.true.
    11073
    … …  
    12992        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
    13093        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
    131                
    132 
    133         write(*,*) "correction for the shape of the ice particles ?"
    134         beta=0.75 ! default value
    135         call getin("ice_shape",beta)
    136         write(*,*) " ice_shape = ",beta
    137        
    138         write(*,*) "time subsampling for microphysic ?"
    139         imicro = 1
    140         call getin("imicro",imicro)
    141         write(*,*)"imicro = ",imicro
    142        
    143         microtimestep = ptimestep/float(imicro)
    144         write(*,*)"Physical timestep is",ptimestep
    145         write(*,*)"Microphysics timestep is",microtimestep
    14694
    14795        firstcall=.false.
    14896      ENDIF ! of IF (firstcall)
    14997     
    150 c-----Initialization
    151       subpdq(1:ngrid,1:nlay,1:nq) = 0
    152       subpdt(1:ngrid,1:nlay)      = 0
    153       pdqscloud(1:ngrid,1:nq)     = 0
    154       zq(1:ngrid,1:nlay,1:nq)     = 0
    15598
    156 c-----Computing the different layer properties for clouds sedimentation     
    157       do l=1,nlay
    158         do ig=1, ngrid
    159           masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
    160           epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
    161          enddo
    162       enddo
    163 
    164 c------------------------------------------------------------------
    165 c Time subsampling for coupled microphysics and sedimentation
    166 c------------------------------------------------------------------
    167       DO microstep=1,imicro
    168      
    169 c-------------------------------------------------------------------
    170 c   1.  Tendencies:
    171 c------------------
    172 
    173 
    174 c------ Temperature tendency subpdt
    175         ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
    176         ! If imicro=1 subpdt is the same as pdt
    177         DO l=1,nlay
    178           DO ig=1,ngrid
    179              subpdt(ig,l) = subpdt(ig,l)
    180      &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
    181           ENDDO
    182         ENDDO
    183 c------ Traceurs tendencies subpdq
    184 c------ At each micro timestep we add pdq in order to have a stepped entry
    185         IF (microphys) THEN
    186           DO l=1,nlay
    187             DO ig=1,ngrid
    188               subpdq(ig,l,igcm_dust_mass) =
    189      &            subpdq(ig,l,igcm_dust_mass)
    190      &          + pdq(ig,l,igcm_dust_mass)
    191               subpdq(ig,l,igcm_dust_number) =
    192      &            subpdq(ig,l,igcm_dust_number)
    193      &          + pdq(ig,l,igcm_dust_number)
    194               subpdq(ig,l,igcm_ccn_mass) =
    195      &            subpdq(ig,l,igcm_ccn_mass)
    196      &          + pdq(ig,l,igcm_ccn_mass)
    197               subpdq(ig,l,igcm_ccn_number) =
    198      &            subpdq(ig,l,igcm_ccn_number)
    199      &          + pdq(ig,l,igcm_ccn_number)
    200             ENDDO
    201           ENDDO
    202         ENDIF
    203         DO l=1,nlay
    204           DO ig=1,ngrid
    205             subpdq(ig,l,igcm_h2o_ice) =
    206      &          subpdq(ig,l,igcm_h2o_ice)
    207      &        + pdq(ig,l,igcm_h2o_ice)
    208             subpdq(ig,l,igcm_h2o_vap) =
    209      &          subpdq(ig,l,igcm_h2o_vap)
    210      &        + pdq(ig,l,igcm_h2o_vap)
    211           ENDDO
    212         ENDDO
    213        
    214        
    215 c-------------------------------------------------------------------
    216 c   2.  Main call to the different cloud schemes:
    217 c------------------------------------------------
    218         IF (microphys) THEN
    219            CALL improvedclouds(ngrid,nlay,microtimestep,
    220      &             pplev,pplay,pzlev,pt,subpdt,
    221      &             pq,subpdq,subpdqcloud,subpdtcloud,
     99c     Main call to the different cloud schemes:
     100      IF (microphys) THEN
     101        CALL improvedclouds(ngrid,nlay,ptimestep,
     102     &             pplev,pplay,pt,pdt,
     103     &             pq,pdq,pdqcloud,pdtcloud,
    222104     &             nq,tauscaling,rdust,rice,nuice,
    223105     &             rsedcloud,rhocloud)
    224         ELSE
    225            CALL simpleclouds(ngrid,nlay,microtimestep,
    226      &             pplev,pplay,pzlev,pzlay,pt,subpdt,
    227      &             pq,subpdq,subpdqcloud,subpdtcloud,
     106      ELSE
     107        CALL simpleclouds(ngrid,nlay,ptimestep,
     108     &             pplev,pplay,pzlev,pzlay,pt,pdt,
     109     &             pq,pdq,pdqcloud,pdtcloud,
    228110     &             nq,tau,rice,nuice,rsedcloud)
    229         ENDIF
    230      
    231 c--------------------------------------------------------------------
    232 c    3.  CCN & ice sedimentation:
    233 c--------------------------------
    234 ! Sedimentation is done here for water ice and its CCN only
    235 ! callsedim in physics is done for dust (and others species if any)           
    236                    
    237         DO l=1,nlay
    238          DO ig=1,ngrid
    239           subpdt(ig,l) =
    240      &      subpdt(ig,l) + subpdtcloud(ig,l)
    241          ENDDO
    242         ENDDO
    243        
    244 c------- water ice update before sedimentation (radius is done in the cloud scheme routine)   
    245          DO l=1,nlay
    246           DO ig=1, ngrid
    247           zq(ig,l,igcm_h2o_ice)= max(1e-30,
    248      &       pq(ig,l,igcm_h2o_ice) + (subpdq(ig,l,igcm_h2o_ice)
    249      &       + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep)
    250 !          zq(ig,l,igcm_h2o_vap)= max(1e-30,
    251 !     &       pq(ig,l,igcm_h2o_vap) + (subpdq(ig,l,igcm_h2o_vap)
    252 !     &       + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep)
    253           ENDDO
    254          ENDDO
    255        
    256      
    257 c------- CCN update before sedimentation 
    258         IF (microphys) THEN
    259          DO l=1,nlay
    260           DO ig=1,ngrid
    261            zq(ig,l,igcm_ccn_number)=
    262      &       pq(ig,l,igcm_ccn_number) + (subpdq(ig,l,igcm_ccn_number)
    263      &       + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep
    264            zq(ig,l,igcm_ccn_number)=  max(1e-30,
    265      &       zq(ig,l,igcm_ccn_number))!*tauscaling(ig)) ! OU pas tauscaling ?
    266            zq(ig,l,igcm_ccn_mass)=
    267      &       pq(ig,l,igcm_ccn_mass) + (subpdq(ig,l,igcm_ccn_mass)
    268      &       + subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep
    269            zq(ig,l,igcm_ccn_mass)=max(1e-30,
    270      &       zq(ig,l,igcm_ccn_mass))!*tauscaling(ig)) ! OU pas tauscaling ?
    271           ENDDO
    272          ENDDO
    273         ENDIF
    274        
    275 
    276        
    277         IF (microphys) THEN
    278        
    279 c------- CCN number sedimentation
    280           call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    281      &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
    282      &        rhocloud,zq(1,1,igcm_ccn_number),wq,beta)
    283           do ig=1,ngrid
    284             ! matters if one would like to know ccn surface deposition
    285             pdqscloud(ig,igcm_ccn_number)=
    286      &          pdqscloud(ig,igcm_ccn_number) + wq(ig,1)
    287           enddo
    288          
    289 c------- CCN mass sedimentation
    290           call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    291      &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
    292      &        rhocloud,zq(1,1,igcm_ccn_mass),wq,beta)
    293           do ig=1,ngrid
    294             ! matters if one would like to know ccn surface deposition
    295             pdqscloud(ig,igcm_ccn_mass)=
    296      &          pdqscloud(ig,igcm_ccn_mass) + wq(ig,1)
    297           enddo
    298          
    299 c------- Water ice sedimentation -- improved microphys
    300           call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    301      &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
    302      &        rhocloud,zq(1,1,igcm_h2o_ice),wq,beta)
    303         ELSE
    304        
    305 c------- Water ice sedimentation -- simple microphys
    306           call newsedim(ngrid,nlay,ngrid*nlay,1,
    307      &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
    308      &        rho_q(igcm_h2o_ice),zq(1,1,igcm_h2o_ice),wq,beta)
    309      
    310         ENDIF
    311 
    312 
    313 c------- Surface ice tendency update
    314         DO ig=1,ngrid
    315           pdqscloud(ig,igcm_h2o_ice)=
    316      &          pdqscloud(ig,igcm_h2o_ice) + wq(ig,1)
    317         ENDDO
    318      
    319 
    320 c-------------------------------------------------------------------
    321 c   5.  Updating tendencies after sedimentation:
    322 c-----------------------------------------------
    323 
    324         DO l=1,nlay
    325          DO ig=1,ngrid
    326            
    327            subpdq(ig,l,igcm_h2o_ice) =
    328      &     (zq(ig,l,igcm_h2o_ice) - pq(ig,l,igcm_h2o_ice))
    329      &      /microtimestep
    330 
    331      
    332            subpdq(ig,l,igcm_h2o_vap)=subpdq(ig,l,igcm_h2o_vap)
    333      &        +subpdqcloud(ig,l,igcm_h2o_vap)
    334      
    335          ENDDO
    336         ENDDO
    337      
    338      
    339         IF (microphys) then
    340          DO l=1,nlay
    341           DO ig=1,ngrid
    342            subpdq(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
    343      &        -pq(ig,l,igcm_ccn_number))/microtimestep
    344            subpdq(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
    345      &        -pq(ig,l,igcm_ccn_mass))/microtimestep
    346           ENDDO
    347          ENDDO
    348         ENDIF
    349        
    350 
    351 
    352       ENDDO ! of DO microstep=1,imicro
    353      
    354 c-------------------------------------------------------------------
    355 c   6.  Compute final tendencies after time loop:
    356 c------------------------------------------------
    357 
    358 c------ Whole temperature tendency pdtcloud
    359        DO l=1,nlay
    360          DO ig=1,ngrid
    361              pdtcloud(ig,l) =
    362      &         subpdt(ig,l)/imicro-pdt(ig,l)
    363           ENDDO
    364        ENDDO
    365 c------ Traceurs tendencies pdqcloud
    366        DO iq=1,nq
    367         DO l=1,nlay
    368          DO ig=1,ngrid
    369             pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/imicro
    370      &       - pdq(ig,l,iq)
    371          ENDDO
    372         ENDDO
    373        ENDDO
    374 c------ Traceurs surface tendencies pdqscloud
    375        DO iq=1,nq
    376         DO ig=1,ngrid
    377            pdqscloud(ig,iq) =
    378      &         pdqscloud(ig,iq)/ptimestep
    379         ENDDO
    380        ENDDO
    381        
    382 
    383          
    384 c------Update the ice particle size "rice" for output or photochemistry.
    385 c------It is not used again in the water cycle.
    386        IF(scavenging) THEN
    387         DO l=1, nlay
    388          DO ig=1,ngrid
    389           Mo = pq(ig,l,igcm_h2o_ice)
    390      &         + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
    391      &         + (pq(ig,l,igcm_ccn_mass)
    392      &          + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
    393      &         *tauscaling(ig) + 1.e-30
    394           No = (pq(ig,l,igcm_ccn_number)
    395      &          + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
    396      &         *tauscaling(ig) + 1.e-30
    397           rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) + 
    398      &                   pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
    399      &                   / Mo * rho_ice
    400      &                  + (pq(ig,l,igcm_ccn_mass)
    401      &                  + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
    402      &                   *tauscaling(ig)/ Mo * rho_dust
    403           rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
    404           if ((Mo.lt.1.e-15) .or. (No.le.50)) then
    405               rice(ig,l) = 1.e-8
    406           else
    407               !! AS: only perform computations if conditions not met
    408               rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
    409           endif
    410          ENDDO
    411         ENDDO
    412        ELSE
    413         DO l=1,nlay
    414           DO ig=1,ngrid
    415             ccntyp =
    416      &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
    417             ccntyp = ccntyp /ccn_factor
    418             rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
    419      &       + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
    420      &      +ccntyp*(4./3.)*pi*rdust(ig,l)*rdust(ig,l)*rdust(ig,l))
    421      &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
    422           ENDDO
    423         ENDDO
    424        ENDIF ! of IF(scavenging)
    425                
    426 !--------------------------------------------------------------
    427 !--------------------------------------------------------------
     111      ENDIF
    428112     
    429113
    … …  
    439123      end do
    440124
    441 c=======================================================================
    442 
     125      RETURN
    443126      END
    444127 
Note: See TracChangeset for help on using the changeset viewer.