Ignore:
Timestamp:
Nov 7, 2011, 6:39:24 PM (13 years ago)
Author:
aslmd
Message:

LMDZ.MARS:

A FIRST VERSION WITH SCAVENGING
TRANSPARENT TO CASUAL USER
FOR WATER CYCLE FOLKS, SEE BELOW

[TN and AS solved conflicts due to r330 >> r357]

07/11/11 == JBM

Changed watercloud.F to call two separate routines,

simpleclouds.F or improvedclouds.F, which are a simplified and
full microphysical scheme for cloud formation, respectively.
Removed the tag called "improved" in watercloud.F, and added
another tag called "microphys" which is defined in callphys.F
instead. Changed routines: callkeys, inifis, physiq, watercloud.

Reimplemented the use of typical profiles of dust particle sizes

and CCNs in simpleclouds.F. Corrected the previously used
analytical CCN profile. Moved ccn_factor to simpleclouds.F,
which won't be used in the new microphysical scheme. Changed
routines: aeropacity, physiq, simpleclouds, watercloud.

Computed rdust at the end of callsedim instead of updatereffrad,

except at firstcall. Renamed rfall into rsedcloud and moved it
in simpleclouds. Moved nuice_sed to initracer.F and added it to
tracer.h. Changed routines: callsedim, physiq, tracer.h,
watercloud, initracer, simpleclouds, updatereffrad.

Added two tracers for the CCNs. Added the different tests in

initracer.F to make sure that, for example, microphys is not
used without doubleq, etc. Corrected an inconsistency in
callsedim.F, and changed the way r0 is updated. Changes
routines: callsedim, inifis, initracer, physiq, testphys1d,
tracer.h.

Added the ability to have a spatially variable density in

newsedim (same method as that used for the radius of
sedimentation). Required the addition of one input to newsedim,
which is the size of the density variable "rho". Changed
routines: callsedim, newsedim.

Added an output to aeropacity called "tauscaling", which is a

factor to convert dust_mass and dust_number to true quantities,
based on tauvis. Removed ccn and qdust from aeropacity, which
are obsolete now.

Wrote improvedcloud.F which includes all the microphysical

scheme, and connected it to the sedimentation scheme. Added and
changed routines: callsedim, physiq, growthrate, nuclea,
improvedclouds, initracer, watercloud, microphys.h.

07/11/11 == TN

Added CCN & dust tracers updates in physiq.F

Corrected a bug that can give negative CCN numbers, due to the
use of single precision values (only 6 decimals) whereas up to 10E+10
CCN can be formed or disappears...

Corrected physical bug that causes h2o_vap or h2o_ice

to be negative in improvedclouds.F.

Corrected physical bug that causes CCN & dust tracers

to be negative in improvedclouds.F when all ice sublimates and
releases dust

Added parameter contact mteta in callphys.def

Default value is still 0.95, see inifis.F

Changed tendancy computation for dust_number in improvedclouds.F

that was not the right one. Indeed, the scavenged dust_number tracer
is computed from the dust_mass one, and its tendancy before scavenging
must be taken into account to compute its scavenging's tendancy.

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
5 added
13 edited

Legend:

Unmodified
Added
Removed

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

    r222 r358  
    11      SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
    2      &    pq,ccn,tauref,tau,aerosol,reffrad,nueffrad,
     2     &    pq,tauscaling,tauref,tau,aerosol,reffrad,nueffrad,
    33     &    QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
    44                                                   
    … …  
    9393      REAL taudusttmp(ngridmx)! Temporary dust opacity
    9494                               !   used before scaling
     95      REAL tauscaling(ngridmx) ! Scaling factor for qdust and Ndust
    9596      REAL taudustvis(ngridmx) ! Dust opacity after scaling
    9697      REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as
    … …  
    102103                               !   Qabs instead of Qext
    103104                               !   (direct comparison with TES)
    104       REAL qdust(ngridmx,nlayermx) ! True dust mass mixing ratio
    105       REAL ccn(ngridmx,nlayermx)   ! Cloud condensation nuclei
    106                                    !   (particules kg-1)
    107       REAL qtot(ngridmx)           ! Dust column (kg m-2)
    108 
    109 c     CCN reduction factor
    110       REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant
    111 
    112 c
     105
    113106c   local saved variables
    114107c   ---------------------
    … …  
    194187!       otherwise default value read from starfi.nc file will be used)
    195188        call getin("tauvis",tauvis)
    196 
    197 !       Some information about the water cycle
    198         IF (water) THEN
    199           write(*,*) "water_param CCN reduc. fac. ", ccn_factor
    200         ENDIF
    201189
    202190        firstcall=.false.
    … …  
    424412c -----------------------------------------------------------------
    425413
     414c     Temporary scaling factor
    426415      taudusttmp(1:ngrid)=0.
    427416      DO iaer=1,naerdust
    … …  
    434423        ENDDO
    435424      ENDDO
     425
     426c     Saved scaling factor
     427      DO ig=1,ngrid
     428          tauscaling(ig) = tauref(ig) *
     429     &                     pplev(ig,1) / 700.E0 / taudusttmp(ig)
     430      ENDDO
     431
     432c     Opacity computation
    436433      DO iaer=1,naerdust
    437434        DO l=1,nlayer
    … …  
    446443        ENDDO
    447444      ENDDO
    448 
    449 c -----------------------------------------------------------------
    450 c Computing the number of condensation nuclei
    451 c -----------------------------------------------------------------
    452       DO iaer = 1, naerkind ! Loop on aerosol kind
    453 c     --------------------------------------------
    454         aerkind2: SELECT CASE (name_iaer(iaer))
    455 c==================================================================
    456         CASE("dust_conrath") aerkind2     ! Typical dust profile
    457 c==================================================================
    458           DO l=1,nlayer
    459             DO ig=1,ngrid
    460               ccn(ig,l) = max(aerosol(ig,l,iaer) /
    461      &                  pi / QREFvis3d(ig,l,iaer) *
    462      &                  (1.+nueffrad(ig,l,iaer))**3. /
    463      &                  reffrad(ig,l,iaer)**2. * g /
    464      &                  (pplev(ig,l)-pplev(ig,l+1)),1e-30)
    465             ENDDO
    466           ENDDO
    467 c==================================================================
    468         CASE("dust_doubleq") aerkind2!Two-moment scheme for dust
    469 c        (transport of mass and number mixing ratio)
    470 c==================================================================
    471           qtot(1:ngridmx) = 0.
    472           DO l=1,nlayer
    473             DO ig=1,ngrid
    474               qdust(ig,l) = pq(ig,l,igcm_dust_mass) * tauref(ig) *
    475      &                      pplev(ig,1) / 700.E0 / taudusttmp(ig)
    476               qtot(ig) = qtot(ig) + qdust(ig,l) *
    477      &                   (pplev(ig,l)-pplev(ig,l+1)) / g
    478               ccn(ig,l) = max( ( ref_r0 /
    479      &                    reffrad(ig,l,iaer) )**3. *
    480      &                    r3n_q * qdust(ig,l) ,1e-30)
    481             ENDDO
    482           ENDDO
    483 c==================================================================
    484         END SELECT aerkind2
    485 c     -----------------------------------
    486       ENDDO ! iaer (loop on aerosol kind)
    487 
    488 
    489 c -----------------------------------------------------------------
    490 c -----------------------------------------------------------------
    491 c  Reduce number of nuclei
    492 !         TEMPORAIRE : r�duction du nombre de nuclei FF 04/200
    493 !         reduction facteur 3
    494 !         ccn(ig,l) = ccn(ig,l) / 27.
    495 !         reduction facteur 2
    496 !         ccn(ig,l) = ccn(ig,l) / 8.
    497 c -----------------------------------------------------------------
    498        DO l=1,nlayer
    499          DO ig=1,ngrid
    500             ccn(ig,l) = ccn(ig,l) / ccn_factor
    501          ENDDO
    502        ENDDO
    503 c -----------------------------------------------------------------
    504 c -----------------------------------------------------------------
    505 
    506445
    507446c -----------------------------------------------------------------
    … …  
    548487     &               'm2.kg-1',3,dsodust)
    549488        ENDIF
    550 c       CALL WRITEDIAGFI(ngridmx,'ccn','Cond. nuclei',
    551 c    &                    'part kg-1',3,ccn)
    552489      ELSE
    553         CALL writeg1d(ngrid,nlayer,dsodust,'dsodust','m2.kg-1')
    554       ENDIF
     490        CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
     491     &                   dsodust)
     492      ENDIF ! of IF (ngrid.NE.1)
    555493c -----------------------------------------------------------------
    556494      return

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r284 r358  
    1212     &   ,callg2d,linear,rayleigh,tracer,active,doubleq,submicron       &
    1313     &   ,lifting,callddevil,scavenging,sedimentation,activice,water    &
    14      &   ,caps,photochem,calltherm,outptherm,callrichsl,callslope
     14     &   ,microphys,caps,photochem,calltherm,outptherm,callrichsl       &
     15     &   ,callslope
    1516     
    1617      COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
    … …  
    4849      integer dustbin
    4950      logical active,doubleq,submicron,lifting,callddevil,scavenging
    50       logical sedimentation,activice,water,caps
     51      logical sedimentation,activice,water,microphys,caps
    5152      logical photochem
    5253      integer nqchem_min

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

    r353 r358  
    22     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
    33     $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
     4     &     tauref,tau,aerosol,tauscaling,rdust,rice,nuice,co2ice)
    45     &     tauref,tau,aerosol,ccn,rdust,rice,nuice,co2ice)
    56
    … …  
    160161
    161162      REAL pq(ngrid,nlayer,nq)
    162       REAL ccn(ngridmx,nlayermx)   ! Cloud condensation nuclei
    163                                    !   (particules kg-1)
     163      REAL tauscaling(ngridmx) ! Conversion factor for
     164                               ! qdust and Ndust
    164165      REAL albedo(ngrid,2),emis(ngrid)
    165166      REAL ls,zday
    … …  
    395396c     Computing aerosol optical depth in each layer:
    396397      CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
    397      &      pq,ccn,tauref,tau,aerosol,reffrad,nueffrad,
     398     &      pq,tauscaling,tauref,tau,aerosol,reffrad,nueffrad,
    398399     &      QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
    399400

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

    r82 r358  
    11      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
    2      $                pplev,zlev, pt, rdust, rice,
     2     &                pplev,zlev, pt, rdust, rice,
     3     &                rsedcloud,rhocloud,
    34     &                pq, pdqfi, pdqsed,pdqs_sed,nq)
    45      IMPLICIT NONE
    … …  
    5960      real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
    6061      real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
    61       real r0(ngridmx,nlayermx) ! geometric mean doubleq radius (m)
     62      real r0(ngridmx,nlayermx) ! geometric mean radius used for
     63                                !   sedimentation (m)
     64      real r0dust(ngridmx,nlayermx) ! geometric mean radius used for
     65                                    !   dust (m)
     66      real r0ccn(ngridmx,nlayermx)  ! geometric mean radius used for
     67                                    !   CCNs (m)
    6268c     Sedimentation radius of water ice
    63       real rfall(ngridmx,nlayermx)
    64 c     Sedimentation effective variance of water ice
    65       REAL, PARAMETER :: nuice_sed = 0.45  !! TESTS_JB  !! 0.1 avant
     69      real rsedcloud(ngridmx,nlayermx)
     70c     Cloud density (kg.m-3)
     71      real rhocloud(ngridmx,nlayermx)
     72
    6673
    6774c     Discrete size distributions (doubleq)
    … …  
    8996      save rr
    9097      integer radpower
     98      real sigma0
    9199
    92100c       3) Other local variables used in doubleq
    93101
    94       real reff(ngridmx,nlayermx,2)  ! for diagnostic only
    95102      INTEGER idust_mass   ! index of tracer containing dust mass
    96103                           !   mix. ratio
    97104      INTEGER idust_number ! index of tracer containing dust number
    98105                           !   mix. ratio
     106      INTEGER iccn_mass    ! index of tracer containing CCN mass
     107                           !   mix. ratio
     108      INTEGER iccn_number  ! index of tracer containing CCN number
     109                           !   mix. ratio
    99110      SAVE idust_mass,idust_number
     111      SAVE iccn_mass,iccn_number
    100112
    101113c     Firstcall:
    … …  
    160172        ENDIF !of if (doubleq)
    161173
     174        IF (scavenging) THEN
     175          iccn_mass=0
     176          iccn_number=0
     177          do iq=1,nq
     178            if (noms(iq).eq."ccn_mass") then
     179              iccn_mass=iq
     180            endif
     181            if (noms(iq).eq."ccn_number") then
     182              iccn_number=iq
     183            endif
     184          enddo
     185          ! check that we did find the tracers
     186          if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
     187            write(*,*) 'callsedim: error! could not identify'
     188            write(*,*) ' tracers for ccn mass and number mixing'
     189            write(*,*) ' ratio and scavenging is activated!'
     190            stop
     191          endif
     192        ENDIF !of if (scavenging)
     193
    162194        IF (water) THEN
    163195          write(*,*) "water_param nueff Sedimentation:", nuice_sed
    … …  
    199231
    200232c =================================================================
     233c     Compute the geometric mean radius used for sedimentation
     234
     235      if (doubleq) then
     236        do l=1,nlay
     237          do ig=1, ngrid
     238            r0dust(ig,l) =
     239     &        CBRT(r3n_q*zqi(ig,l,idust_mass)/
     240     &        max(zqi(ig,l,idust_number),0.01))
     241            r0dust(ig,l)=min(max(r0dust(ig,l),1.e-10),500.e-6)
     242          end do
     243        end do
     244      endif
     245      if (scavenging) then
     246        do l=1,nlay
     247          do ig=1, ngrid
     248            r0ccn(ig,l) = rsedcloud(ig,l)/(1.+nuice_sed)**4.5
     249          end do
     250        end do
     251      endif
     252
     253c =================================================================
    201254      do iq=1,nq
    202255        if(radius(iq).gt.1.e-9) then   ! no sedim for gaz
    … …  
    205258c         DOUBLEQ CASE
    206259c -----------------------------------------------------------------
    207           if (doubleq.and.
     260          if ((doubleq.and.
    208261     &        ((iq.eq.idust_mass).or.
    209      &         (iq.eq.idust_number))) then
     262     &         (iq.eq.idust_number))).or.
     263     &        (scavenging.and.
     264     &        ((iq.eq.iccn_mass).or.
     265     &        (iq.eq.iccn_number)))) then
    210266
    211267c           Computing size distribution:
    212268c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    213269
    214             do  l=1,nlay
    215               do ig=1, ngrid
    216                 r0(ig,l)=
    217      &            CBRT(r3n_q*zqi(ig,l,idust_mass)/
    218      &            max(zqi(ig,l,idust_number),0.01))
    219                 r0(ig,l)=min(max(r0(ig,l),1.e-10),500.e-6)
     270            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
     271              do  l=1,nlay
     272                do ig=1, ngrid
     273                  r0(ig,l)=r0dust(ig,l)
     274                end do
    220275              end do
    221             end do
     276              sigma0 = varian
     277            else
     278              do  l=1,nlay
     279                do ig=1, ngrid
     280                  r0(ig,l)=r0ccn(ig,l)
     281                end do
     282              end do
     283              sigma0 = sqrt(log(1.+nuice_sed))
     284            endif
    222285
    223286c        Computing mass mixing ratio for each particle size
    224287c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    225           IF (iq.EQ.idust_mass) then
     288          IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
    226289            radpower = 2
    227290          ELSE
    … …  
    237300                 qr(ig,l,ir)=0.5*(rr(2,ir)-rr(1,ir))*
    238301     &             (rr(1,ir)**radpower)*
    239      &             exp(-(log(rr(1,ir)/r0(ig,l)))**2/(2*varian**2))
     302     &             exp(-(log(rr(1,ir)/r0(ig,l)))**2/(2*sigma0**2))
    240303                 do iint=2,ninter-1
    241304                   qr(ig,l,ir)=qr(ig,l,ir) +
    … …  
    243306     &             (rr(iint,ir)**radpower)*
    244307     &             exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
    245      &             (2*varian**2))
     308     &             (2*sigma0**2))
    246309                 end do
    247310                 qr(ig,l,ir)=qr(ig,l,ir) +
    … …  
    249312     &             (rr(ninter,ir)**radpower)*
    250313     &             exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
    251      &             (2*varian**2))
     314     &             (2*sigma0**2))
    252315
    253316c                **************** old method (not recommended!)
    254317c                qr(ig,l,ir)=(rd(ir)**(5-3*iq))*
    255 c    &           exp( -(log(rd(ir)/r0(ig,l)))**2 / (2*varian**2) )
     318c    &           exp( -(log(rd(ir)/r0(ig,l)))**2 / (2*sigma0**2) )
    256319c                ******************************
    257320
    … …  
    272335c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    273336
    274 c         call zerophys(ngridmx*nlayermx,zqi(1,1,iq))
    275337          zqi(1:ngrid,1:nlay,iq) = 0.
    276 c         call zerophys(ngridmx,pdqs_sed(1,iq))
    277338          pdqs_sed(1:ngrid,iq) = 0.
    278339
    279340          do ir=1,nr
    280              call newsedim(ngrid,nlay,1,ptimestep,
    281      &       pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
    282      &       wq,0.5)
     341             IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then
     342               call newsedim(ngrid,nlay,1,1,ptimestep,
     343     &         pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
     344     &         wq,0.5)
     345             ELSE
     346               call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep,
     347     &         pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir),
     348     &         wq,1.0)
     349             ENDIF
    283350
    284351c            Tendencies
    … …  
    298365c -----------------------------------------------------------------
    299366          else if (water.and.(iq.eq.igcm_h2o_ice)) then
    300 c           if (doubleq) then
    301 c             do l=1,nlay
    302 c               do ig=1,ngrid
    303 c                 rfall(ig,l)=max( rice(ig,l),rdust(ig,l) )
    304 c                 rfall(ig,l)=min(rfall(ig,l),1.e-4)
    305 c               enddo
    306 c             enddo ! of do l=1,nlay
    307 c           else
    308               do l=1,nlay
    309                 do ig=1,ngrid
    310 c               For water cycle, a typical dust size is assumed:
    311 c               r(z)=r0*exp(-z/H) with r0=0.8 micron and H=18 km.
    312 c               rfall(ig,l)=max( rice(ig,l)*1.5,rdust(ig,l) )
    313                 rfall(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3.,
    314      &                           rdust(ig,l) )
    315 c modif FranckMM pour ameliorer cycle H2O: rfall= 20 microns
    316 c mars commente pour l'instant               rfall(ig,l)=20.e-6
    317                 rfall(ig,l)=min(rfall(ig,l),1.e-4)
    318                 enddo
    319               enddo ! of do l=1,nlay
    320 c           endif
    321             call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
    322      &      pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi(1,1,iq),
    323      &      wq,1.0)
     367            if (microphys) then
     368              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,ptimestep,
     369     &        pplev,masse,epaisseur,pt,rsedcloud,rhocloud,zqi(1,1,iq),
     370     &        wq,1.0)
     371            else
     372              call newsedim(ngrid,nlay,ngrid*nlay,1,ptimestep,
     373     &        pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),zqi(1,1,iq),
     374     &        wq,1.0)
     375            endif ! of if (microphys)
    324376c           Tendencies
    325377c           ~~~~~~~~~~
    … …  
    331383c -----------------------------------------------------------------
    332384          else
    333             call newsedim(ngrid,nlay,1,ptimestep,
     385            call newsedim(ngrid,nlay,1,1,ptimestep,
    334386     &      pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
    335387     &      zqi(1,1,iq),wq,1.0)
    … …  
    342394c -----------------------------------------------------------------
    343395
     396c         Compute the final tendency:
     397c         ---------------------------
    344398          DO l = 1, nlay
    345399            DO ig=1,ngrid
    … …  
    353407      enddo ! of do iq=1,nq
    354408 
     409c     Update the dust particle size "rdust"
     410c     -------------------------------------
     411      DO l = 1, nlay
     412        DO ig=1,ngrid
     413          rdust(ig,l)=
     414     &      CBRT(r3n_q*zqi(ig,l,idust_mass)/
     415     &      max(zqi(ig,l,idust_number),0.01))
     416          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
     417        ENDDO
     418      ENDDO
     419
    355420      RETURN
    356421      END

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

    r120 r358  
    1       subroutine growthrate(timestep,t,p,ph2o,psat,seq,r,Cste)
     1      subroutine growthrate(timestep,temp,pmid,ph2o,psat,
     2     &                      seq,rcrystal,growth)
    23
    34      IMPLICIT NONE
    … …  
    78c     Determination of the water ice crystal growth rate
    89c
     10c     Authors: F. Montmessin
     11c       Adapted for the LMD/GCM by J.-B. Madeleine (October 2011)
     12c     
    913c=======================================================================
    1014
    … …  
    1317c   -------------
    1418
     19#include "dimensions.h"
     20#include "dimphys.h"
     21#include "comcstfi.h"
     22#include "tracer.h"
     23#include "microphys.h"
     24
    1525c
    1626c   arguments:
    … …  
    1828
    1929      REAL timestep
    20       REAL t    ! temperature in the middle of the layer (K)
    21       REAL p    ! pressure in the middle of the layer (K)
    22       REAL*8 ph2o ! water vapor partial pressure (Pa)
    23       REAL*8 psat ! water vapor saturation pressure (Pa)
    24       REAL r    ! crystal radius before condensation (m)
    25       REAL*8 seq  ! Equilibrium saturation ratio
    26 !      REAL dr   ! crystal radius variation (m)
    27       REAL*8 Cste
     30      REAL temp     ! temperature in the middle of the layer (K)
     31      REAL pmid     ! pressure in the middle of the layer (K)
     32      REAL*8 ph2o   ! water vapor partial pressure (Pa)
     33      REAL*8 psat   ! water vapor saturation pressure (Pa)
     34      REAL rcrystal ! crystal radius before condensation (m)
     35      REAL*8 seq    ! Equilibrium saturation ratio
     36      REAL*8 growth
    2837
    2938c   local:
    3039c   ------
    3140
    32       REAL*8 molco2,molh2o
    33       REAL*8 Mco2,Mh2o,rho_i,sigh2o
    34       REAL*8 nav,rgp,kbz,pi,To
    35 
    36 c     Effective gas molecular radius (m)
    37       data molco2/2.2d-10/   ! CO2
    38 c     Effective gas molecular radius (m)
    39       data molh2o/1.2d-10/   ! H2O
    40 c     Molecular weight of CO2
    41       data Mco2/44.d-3/            ! kg.mol-1
    42 c     Molecular weight of H2O
    43       data Mh2o/18.d-3/            ! kg.mol-1
    44 c     surface tension of ice/vapor
    45       data sigh2o/0.12/      ! N.m
    46 c     Ice density
    47       data rho_i/917./        ! kg.m-3 also defined in initcld.f
    48 c     Avogadro number
    49       data nav/6.023d23/
    50 c     Perfect gas constant
    51       data rgp/8.3143/
    52 c     Boltzman constant
    53       data kbz/1.381d-23/
    54 c     pi number
    55       data pi/3.141592654/
    56 c     Reference temperature, T=273,15 K
    57       data To/273.15/
    58 
    5941      REAL*8 k,Lv                 
    6042      REAL*8 knudsen           ! Knudsen number (gas mean free path/particle radius)
    61       REAL*8 a,Dv,lambda       ! Intermediate computations for growth rate
     43      REAL*8 afactor,Dv,lambda       ! Intermediate computations for growth rate
    6244      REAL*8 Rk,Rd
    6345
    … …  
    7052
    7153c     - Equilibrium saturation accounting for KeLvin Effect
    72        seq=exp(2*sigh2o*Mh2o/(rho_i*rgp*t*r))
     54c      seq=exp(2*sigh2o*mh2o/(rho_ice*rgp*t*r))
     55c      (already computed in improvedcloud.F)
    7356
    7457c     - Thermal conductibility of CO2
    75       k  = (0.17913 * t - 13.9789) * 4.184e-4
     58      k  = (0.17913 * temp - 13.9789) * 4.184e-4
    7659c     - Latent heat of h2o (J.kg-1)
    77       Lv = (2834.3 - 0.28 * (t-To) - 0.004 * (t-To)**2 ) * 1.e+3
     60      Lv = (2834.3 - 0.28 * (temp-To) - 0.004 * (temp-To)**2 ) * 1.e+3
    7861
    7962c     - Constant to compute gas mean free path
    8063c     l= (T/P)*a, with a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
    81       a = 0.707*rgp/(4 * pi* molco2**2  * nav)
     64      afactor = 0.707*rgp/(4 * pi* molco2**2  * nav)
    8265
    8366c     - Compute Dv, water vapor diffusion coefficient
    … …  
    8568c       the nature of which depending on the Knudsen number.
    8669
    87       Dv = 1./3. * sqrt( 8*kbz*t/(pi*Mh2o/nav) )* kbz * t /
    88      &   ( pi * p * (molco2+molh2o)**2 * sqrt(1.+Mh2o/Mco2) )
     70      Dv = 1./3. * sqrt( 8*kbz*temp/(pi*mh2o/nav) )* kbz * temp /
     71     &   ( pi * pmid * (molco2+molh2o)**2 * sqrt(1.+mh2o/mco2) )
    8972
    90       knudsen = t / p * a / r
     73      knudsen = temp / pmid * afactor / rcrystal
    9174      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen)
    9275      Dv      = Dv / (1. + lambda * knudsen)
    9376
    9477c     - Compute Rk
    95       Rk = Lv**2 * rho_i * Mh2o / (k*rgp*t**2.)
     78      Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp**2.)
    9679c     - Compute Rd
    97       Rd = rgp * t *rho_i / (Dv*psat*Mh2o)
     80      Rd = rgp * temp *rho_ice / (Dv*psat*mh2o)
    9881
    99 c     - Compute Cste=rdr/dt, then r(t+1)= sqrt(r(t)**2.+2.*Cste*dt)
    100        Cste = 1. / (seq*Rk+Rd)
    101 c      Cste = (ph2o/psat-seq) / (seq*Rk+Rd)
    102 c      rf   = sqrt( max( r**2.+2.*Cste*timestep , 0. ) )
     82c     - Compute growth=rdr/dt, then r(t+1)= sqrt(r(t)**2.+2.*growth*dt)
     83      growth = 1. / (seq*Rk+Rd)
     84c       growth = (ph2o/psat-seq) / (seq*Rk+Rd)
     85c      rf   = sqrt( max( r**2.+2.*growth*timestep , 0. ) )
    10386c      dr   = rf-r
    10487

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

    r288 r358  
    6161#include "datafile.h"
    6262#include "slope.h"
     63#include "microphys.h"
    6364#ifdef MESOSCALE
    6465#include "comsoil.h"     !!MESOSCALE -- needed to fill volcapa
    … …  
    304305! TRACERS:
    305306
     307! dustbin
    306308         write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
    307309         dustbin=0 ! default value
    308310         call getin("dustbin",dustbin)
    309311         write(*,*)" dustbin = ",dustbin
    310 
     312! active
    311313         write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
    312314         active=.false. ! default value
    … …  
    321323           stop
    322324         endif
    323 
     325! doubleq
    324326         write(*,*)"use mass and number mixing ratios to predict",
    325327     &             " dust size ?"
    … …  
    327329         call getin("doubleq",doubleq)
    328330         write(*,*)" doubleq = ",doubleq
    329 
     331! submicron
    330332         submicron=.false. ! default value
    331333         call getin("submicron",submicron)
    … …  
    345347           stop
    346348         endif
    347 
     349! lifting
    348350         write(*,*)"dust lifted by GCM surface winds ?"
    349351         lifting=.false. ! default value
    … …  
    358360           stop
    359361         endif
    360 
     362! callddevil
    361363         write(*,*)" dust lifted by dust devils ?"
    362364         callddevil=.false. !default value
    363365         call getin("callddevil",callddevil)
    364366         write(*,*)" callddevil = ",callddevil
    365          
    366367
    367368! Test of incompatibility:
    … …  
    372373           stop
    373374         endif
    374 
    375          write(*,*)"Dust scavenging by CO2 snowfall ?"
    376          scavenging=.false. ! default value
    377          call getin("scavenging",scavenging)
    378          write(*,*)" scavenging = ",scavenging
    379          
    380 
    381 ! Test of incompatibility:
    382 ! if scavenging is used, then dustbin should be > 0
    383 
    384          if (scavenging.and.(dustbin.lt.1)) then
    385            print*,'if scavenging is used, then dustbin should > 0'
    386            stop
    387          endif
    388 
     375! sedimentation
    389376         write(*,*) "Gravitationnal sedimentation ?"
    390377         sedimentation=.true. ! default value
    391378         call getin("sedimentation",sedimentation)
    392379         write(*,*) " sedimentation = ",sedimentation
    393 
     380! activice
    394381         write(*,*) "Radiatively active transported atmospheric ",
    395382     &              "water ice ?"
    … …  
    397384         call getin("activice",activice)
    398385         write(*,*) " activice = ",activice
    399 
     386! water
    400387         write(*,*) "Compute water cycle ?"
    401388         water=.false. ! default value
    … …  
    412399         if (water.and..not.tracer) then
    413400           print*,'if water is used, tracer should be used too'
     401           stop
     402         endif
     403! microphys
     404         write(*,*)"Microphysical scheme for water-ice clouds?"
     405         microphys=.false. ! default value
     406         call getin("microphys",microphys)
     407         write(*,*)" microphys = ",microphys
     408
     409! microphysical parameter contact       
     410         write(*,*) "water contact parameter ?"
     411         mteta  = 0.95
     412         call getin("mteta",mteta)
     413         write(*,*) "mteta = ", mteta
     414
     415! scavenging
     416         write(*,*)"Dust scavenging by H2O/CO2 snowfall ?"
     417         scavenging=.false. ! default value
     418         call getin("scavenging",scavenging)
     419         write(*,*)" scavenging = ",scavenging
     420         
     421
     422! Test of incompatibility:
     423! if scavenging is used, then dustbin should be > 0
     424
     425         if (scavenging.and.(dustbin.lt.1)) then
     426           print*,'if scavenging is used, then dustbin should > 0'
     427           stop
     428         endif
     429         if ((microphys.and..not.doubleq).or.
     430     &       (microphys.and..not.active).or.
     431     &       (microphys.and..not.scavenging).or.
     432     &       (microphys.and..not.water)) then
     433           print*,'if microphys is used, then doubleq,'
     434           print*,'active, water and scavenging must all be used!'
     435           stop
     436         endif
     437         if ((scavenging.and..not.doubleq).or.
     438     &       (scavenging.and..not.active).or.
     439     &       (scavenging.and..not.water).or.
     440     &       (scavenging.and..not.microphys)) then
     441           print*,'if scavenging is used, then doubleq,'
     442           print*,'active, water and microphys must all be used!'
    414443           stop
    415444         endif
    … …  
    425454         write(*,*) " caps = ",caps
    426455
    427 
     456! albedo_h2o_ice
    428457         write(*,*) "water ice albedo ?"
    429458         albedo_h2o_ice=0.45
    430459         call getin("albedo_h2o_ice",albedo_h2o_ice)
    431460         write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
    432          
     461! inert_h2o_ice
    433462         write(*,*) "water ice thermal inertia ?"
    434463         inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
    435464         call getin("inert_h2o_ice",inert_h2o_ice)
    436465         write(*,*) " inert_h2o_ice = ",inert_h2o_ice
    437          
     466! frost_albedo_threshold
    438467         write(*,*) "frost thickness threshold for albedo ?"
    439468         frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)

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

    r324 r358  
    213213      igcm_dust_mass=0
    214214      igcm_dust_number=0
     215      igcm_ccn_mass=0
     216      igcm_ccn_number=0
    215217      igcm_dust_submicron=0
    216218      igcm_h2o_vap=0
    … …  
    271273        enddo
    272274      endif ! of if (doubleq)
     275      if (scavenging) then
     276        do iq=1,nqmx
     277          if (noms(iq).eq."ccn_mass") then
     278            igcm_ccn_mass=iq
     279            count=count+1
     280          endif
     281          if (noms(iq).eq."ccn_number") then
     282            igcm_ccn_number=iq
     283            count=count+1
     284          endif
     285        enddo
     286      endif ! of if (scavenging)
    273287      if (submicron) then
    274288        do iq=1,nqmx
    … …  
    481495      rho_ice=920.    ! Water ice density (kg.m-3)
    482496      nuice_ref=0.1   ! Effective variance nueff of the
    483                       ! water-ice size distributions
     497                      ! water-ice size distribution
     498      nuice_sed=0.45   ! Sedimentation effective variance
     499                      ! of the water-ice size distribution
    484500
    485501      if (doubleq) then
    … …  
    530546        write(*,*) "initracer: doubleq_param alpha_lift:",
    531547     &    alpha_lift(igcm_dust_mass)
    532 
    533548      else
    534549
    … …  
    549564       endif
    550565      end if    ! (doubleq)
     566
     567
     568c     Scavenging of dust particles by H2O clouds:
     569c     ------------------------------------------
     570c     Initialize the two tracers used for the CCNs
     571      if (water.AND.doubleq.AND.scavenging) then
     572        radius(igcm_ccn_mass) = radius(igcm_dust_mass)
     573        alpha_lift(igcm_ccn_mass) = 1e-30
     574        alpha_devil(igcm_ccn_mass) = 1e-30
     575        rho_q(igcm_ccn_mass) = rho_dust
     576
     577        radius(igcm_ccn_number) = radius(igcm_ccn_mass)
     578        alpha_lift(igcm_ccn_number) = alpha_lift(igcm_ccn_mass)
     579        alpha_devil(igcm_ccn_number) = alpha_devil(igcm_ccn_mass)
     580        rho_q(igcm_ccn_number) = rho_q(igcm_ccn_mass)
     581      endif ! of if (water.AND.doubleq.AND.scavenging)
    551582
    552583c     Submicron dust mode:
    … …  
    670701       endif
    671702
     703       if (scavenging) then
     704       ! verify that we indeed have ccn_mass and ccn_number tracers
     705         if (igcm_ccn_mass.eq.0) then
     706           write(*,*) "initracer: error !!"
     707           write(*,*) "  cannot use scavenging option without ",
     708     &                "a ccn_mass tracer !"
     709           stop
     710         endif
     711         if (igcm_ccn_number.eq.0) then
     712           write(*,*) "initracer: error !!"
     713           write(*,*) "  cannot use scavenging option without ",
     714     &                "a ccn_number tracer !"
     715           stop
     716         endif
     717       endif ! of if (scavenging)
     718
    672719       if (photochem .or. callthermos) then
    673720       ! verify that we indeed have the chemistry tracers

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

    r147 r358  
    1       SUBROUTINE newsedim(ngrid,nlay,naersize,ptimestep,
     1      SUBROUTINE newsedim(ngrid,nlay,naersize,nrhosize,ptimestep,
    22     &  pplev,masse,epaisseur,pt,rd,rho,pqi,wq,beta)
    33      IMPLICIT NONE
    … …  
    2121c   ----------
    2222
    23       INTEGER,INTENT(IN) :: ngrid,nlay,naersize
     23      INTEGER,INTENT(IN) :: ngrid,nlay,naersize,nrhosize
    2424      REAL,INTENT(IN) :: ptimestep            ! pas de temps physique (s)
    2525      REAL,INTENT(IN) :: pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
    … …  
    2828      real,intent(in) :: epaisseur (ngrid,nlay)  ! epaisseur d'une couche (m)
    2929      real,intent(in) :: rd(naersize)             ! particle radius (m)
    30       real,intent(in) :: rho             ! particle density (kg.m-3)
     30      real,intent(in) :: rho(nrhosize)             ! particle density (kg.m-3)
    3131
    3232
    … …  
    4444
    4545      INTEGER l,ig, k, i
    46       REAL rfall
     46      REAL rfall,rhofall
    4747
    4848      LOGICAL,SAVE :: firstcall=.true.
    … …  
    109109        do  l=1,nlay
    110110          do ig=1, ngrid
     111c           radius
    111112            if (naersize.eq.1) then
    112113              rfall=rd(1)
    … …  
    115116              rfall=rd(i)
    116117            endif 
    117             vstokes(ig,l) = b * rho * rfall**2 *
     118c           density
     119            if (nrhosize.eq.1) then
     120              rhofall=rho(1)
     121            else
     122              i=ngrid*(l-1)+ig
     123              rhofall=rho(i)
     124            endif 
     125c           vstokes
     126            vstokes(ig,l) = b * rhofall * rfall**2 *
    118127     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
    119128

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

    r353 r358  
    198198      REAL nuice(ngridmx,nlayermx)   ! Estimated effective variance
    199199                                     !   of the size distribution
     200      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
     201      real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
    200202      REAL surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3, if photochemistry)
    201203      REAL surfice(ngridmx,nlayermx)  !  ice surface area (micron2/cm3, if photochemistry)
    … …  
    231233      REAL fluxtop_sw(ngridmx,2)     !Outgoing SW (solar) flux to space (W.m-2)
    232234      REAL tauref(ngridmx)           ! Reference column optical depth at 700 Pa
    233                                      ! (used if active=F)
    234235      REAL tau(ngridmx,naerkind)     ! Column dust optical depth at each point
    235236      REAL zls                       !  solar longitude (rad)
    … …  
    294295      character*5 str5
    295296      real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
    296       REAL ccn(ngridmx,nlayermx)   ! Cloud condensation nuclei
    297                                    !   (particules kg-1)
    298       SAVE ccn  !! in case iradia != 1
     297      REAL tauscaling(ngridmx)   ! Convertion factor for qdust and Ndust
     298      SAVE tauscaling            ! in case iradia NE 1
    299299      real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
    300       real qtot1,qtot2 ! total aerosol mass
    301300      integer igmin, lmin
    302301      logical tdiag
    … …  
    316315      REAL Qabsice                ! Water ice absorption coefficient
    317316
     317c Test 1d scavenging
     318      real h2o_tot
     319      real ccndust_mass(nlayermx)
     320      real ccndust_number(nlayermx)
    318321
    319322      REAL time_phys
    … …  
    548551     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
    549552     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
    550      &     tauref,tau,aerosol,ccn,rdust,rice,nuice,co2ice)
     553     &     tauref,tau,aerosol,tauscaling,rdust,rice,nuice,co2ice)
    551554
    552555c          Outputs for basic check (middle of domain)
    … …  
    920923     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
    921924     &                pq,pdq,zdqcloud,zdqscloud,zdtcloud,
    922      &                nq,naerkind,tau,
    923      &                ccn,rdust,rice,nuice,
    924      &                surfdust,surfice)
    925 
     925     &                nq,tau,tauscaling,rdust,rice,nuice,
     926     &                rsedcloud,rhocloud)
    926927           if (activice) then
    927928c Temperature variation due to latent heat release
    … …  
    935936! increment water vapour and ice atmospheric tracers tendencies
    936937           IF (water) THEN
    937              DO l=1,nlayer
     938          DO l=1,nlayer
    938939               DO ig=1,ngrid
    939                  pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
    940      &                                   zdqcloud(ig,l,igcm_h2o_vap)
    941                  pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
    942      &                                   zdqcloud(ig,l,igcm_h2o_ice)
     940                 pdq(ig,l,igcm_h2o_vap)=
     941     &                     pdq(ig,l,igcm_h2o_vap)+
     942     &                    zdqcloud(ig,l,igcm_h2o_vap)
     943                 pdq(ig,l,igcm_h2o_ice)=
     944     &                    pdq(ig,l,igcm_h2o_ice)+
     945     &                    zdqcloud(ig,l,igcm_h2o_ice)
     946                 IF (scavenging) THEN
     947                   pdq(ig,l,igcm_ccn_mass)=
     948     &                      pdq(ig,l,igcm_ccn_mass)+
     949     &                      zdqcloud(ig,l,igcm_ccn_mass)
     950                   pdq(ig,l,igcm_ccn_number)=
     951     &                      pdq(ig,l,igcm_ccn_number)+
     952     &                      zdqcloud(ig,l,igcm_ccn_number)
     953!!!!!!!!!!!!!!!!!!!!! We need to check that we have Nccn > 0
     954!!!!!!!!!!!!!!!!!!!!! This is due to single preicions roiunding problems     
     955                   if (((pq(ig,l,igcm_ccn_number) +
     956     &                 pdq(ig,l,igcm_ccn_number)*ptimestep)) .le. 0)
     957     &             then
     958                       pdq(ig,l,igcm_ccn_number) =
     959     &                 - pq(ig,l,igcm_ccn_number)/ptimestep
     960                   endif
     961!!!!!!!!!!!!!!!!!!!!!
     962!!!!!!!!!!!!!!!!!!!!!
     963                   pdq(ig,l,igcm_dust_mass)=
     964     &                      pdq(ig,l,igcm_dust_mass)+
     965     &                      zdqcloud(ig,l,igcm_dust_mass)
     966                   pdq(ig,l,igcm_dust_number)=
     967     &                      pdq(ig,l,igcm_dust_number)+
     968     &                      zdqcloud(ig,l,igcm_dust_number)
     969                 ENDIF
    943970               ENDDO
    944971             ENDDO
    … …  
    960987c        --------------
    961988         IF (photochem .or. thermochem) then
     989            PRINT*, 'SURFDUST,SURFICE TO BE IMPLEMENTED. YAAAAAARG.'
    962990            call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
    963991     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
    … …  
    10401068           call callsedim(ngrid,nlayer, ptimestep,
    10411069     &                pplev,zzlev, pt, rdust, rice,
     1070     &                rsedcloud,rhocloud,
    10421071     &                pq, pdq, zdqsed, zdqssed,nq)
    10431072
    … …  
    13781407     &                    "surface h2o_ice","kg/m2",
    13791408     &                    2,qsurf(1,igcm_h2o_ice))
    1380 
     1409c               call wstats(ngrid,'albedo',
     1410c     &                         'albedo',
     1411c     &                         '',2,albedo(1:ngridmx,1))
    13811412               call wstats(ngrid,"mtot",
    13821413     &                    "total mass of water vapor","kg/m2",
    … …  
    15881619     &                       'total mass of water ice',
    15891620     &                       'kg/m2',2,icetot)
    1590 c            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
    1591 c    &                *mugaz/mmol(igcm_h2o_ice)
    1592 c            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
    1593 c    &                       'mol/mol',3,vmr)
     1621            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
     1622     &                *mugaz/mmol(igcm_h2o_ice)
     1623            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
     1624     &                       'mol/mol',3,vmr)
     1625            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
     1626     &                *mugaz/mmol(igcm_h2o_vap)
     1627            call WRITEDIAGFI(ngridmx,'vmr_h2ovap','h2o vap vmr',
     1628     &                       'mol/mol',3,vmr)
    15941629             CALL WRITEDIAGFI(ngridmx,'reffice',
    15951630     &                       'Mean reff',
    15961631     &                       'm',2,rave)
    1597 c            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
    1598 c    &                       'm',3,rice)
     1632            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
     1633     &                       'm',3,rice)
    15991634c            If activice is true, tauTES is computed in aeropacity.F;
    16001635             if (.not.activice) then
    … …  
    16111646                if (watercaptag(ig)) watercapflag(ig) = 1
    16121647             enddo
    1613              CALL WRITEDIAGFI(ngridmx,'watercaptag',
    1614      &                         'Ice water caps',
    1615      &                         '',2,watercapflag)
     1648c            CALL WRITEDIAGFI(ngridmx,'watercaptag',
     1649c    &                         'Ice water caps',
     1650c    &                         '',2,watercapflag)
    16161651            endif
    1617             CALL WRITEDIAGFI(ngridmx,'albedo',
    1618      &                         'albedo',
    1619      &                         '',2,albedo(1:ngridmx,1))
     1652c           CALL WRITEDIAGFI(ngridmx,'albedo',
     1653c    &                         'albedo',
     1654c    &                         '',2,albedo(1:ngridmx,1))
    16201655           endif !(water)
    16211656
    … …  
    16471682c          call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
    16481683           if (doubleq) then
    1649 c            call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
    1650 c    &                       'kg.m-2',2,qsurf(1,1))
    1651 c            call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
    1652 c    &                       'N.m-2',2,qsurf(1,2))
     1684            call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
     1685     &                       'kg.m-2',2,qsurf(1,1))
     1686            call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
     1687     &                       'N.m-2',2,qsurf(1,2))
    16531688c            call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
    16541689c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
    1655 c            call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
    1656 c    &                       'kg.m-2.s-1',2,zdqssed(1,1))
     1690            call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
     1691     &                       'kg.m-2.s-1',2,zdqssed(1,1))
    16571692             call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
    16581693     &                        'm',3,rdust*ref_r0)
    16591694             call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
    16601695     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
    1661 c            call WRITEDIAGFI(ngridmx,'dustN','Dust number',
    1662 c    &                        'part/kg',3,pq(1,1,igcm_dust_number))
     1696             call WRITEDIAGFI(ngridmx,'dustN','Dust number',
     1697     &                        'part/kg',3,pq(1,1,igcm_dust_number))
    16631698#ifdef MESOINI
    16641699             call WRITEDIAGFI(ngridmx,'dustN','Dust number',
    … …  
    16741709             end do
    16751710           endif ! (doubleq)
     1711
     1712           if (scavenging) then
     1713             call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr',
     1714     &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
     1715             call WRITEDIAGFI(ngridmx,'ccnN','CCN number',
     1716     &                        'part/kg',3,pq(1,1,igcm_ccn_number))
     1717           endif ! (scavenging)
     1718
    16761719c          if (submicron) then
    16771720c            call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
    … …  
    18681911         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
    18691912         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
     1913         call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
    18701914! or output in diagfi.nc (for testphys1d)
    18711915         call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
    18721916         call WRITEDIAGFI(ngridmx,'temp','Temperature',
    18731917     &                       'K',1,zt)
    1874 
     1918     
    18751919         if(tracer) then
    18761920c           CALL writeg1d(ngrid,1,tau,'tau','SI')
    … …  
    18801924     &              trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
    18811925            end do
     1926           if (doubleq) then
     1927             call WRITEDIAGFI(ngridmx,'rdust','rdust',
     1928     &                        'm',1,rdust)
     1929           endif
    18821930         end if
     1931         
     1932cccccccccccccccccccccc scavenging outputs 1D TN ccccccccccccccccccc
     1933ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1934
     1935             CALL WRITEDIAGFI(ngridmx,'tauTESap',
     1936     &                         'tau abs 825 cm-1',
     1937     &                         '',1,tauTES)
     1938     
     1939           h2o_tot = qsurf(1,igcm_h2o_ice)
     1940           do l=1,nlayer
     1941             h2o_tot = h2o_tot +
     1942     &           (zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice))
     1943     &                     * (pplev(ig,l) - pplev(ig,l+1)) / g
     1944             ccndust_mass(l) =
     1945     &       pq(1,l,igcm_dust_mass)+pq(1,l,igcm_ccn_mass)
     1946             ccndust_number(l) =
     1947     &       pq(1,l,igcm_dust_number)+pq(1,l,igcm_ccn_number)
     1948     
     1949           end do
     1950         
     1951             call WRITEDIAGFI(ngrid,'ccn+dust_mass','CCN+Dust mass',
     1952     &                        'kg/kg',1,ccndust_mass)
     1953             call WRITEDIAGFI(ngrid,'ccn+dust_number',
     1954     &       'CCN+Dust number','part/kg',1,ccndust_number)
     1955             call WRITEDIAGFI(ngridmx,'h2o_ice_s',
     1956     &                       'surface h2o_ice',
     1957     &                       'kg.m-2',0,qsurf(1,igcm_h2o_ice))
     1958             call WRITEDIAGFI(ngrid,'h2otot',
     1959     &       'total water mass','kg.m-2',0,h2o_tot)
     1960     
     1961ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1962ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1963
    18831964
    18841965         zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g

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

    r283 r358  
    312312            close(91)
    313313          endif ! of if (txt.eq."dust_number")
     314          ! NB: some more initializations (chemistry) is done later
     315          ! CCN MASS
     316          if (txt.eq."ccn_mass") then
     317            !look for a "profile_ccn_mass" input file
     318            open(91,file='profile_ccn_mass',status='old',
     319     &       form='formatted',iostat=ierr)
     320            if (ierr.eq.0) then
     321              read(91,*) qsurf(iq)
     322              do ilayer=1,nlayermx
     323                read(91,*) q(ilayer,iq)
     324              enddo
     325            else
     326              write(*,*) "No profile_ccn_mass file!"
     327            endif
     328            close(91)
     329          endif ! of if (txt.eq."ccn_mass")
     330          ! CCN NUMBER
     331          if (txt.eq."ccn_number") then
     332            !look for a "profile_ccn_number" input file
     333            open(91,file='profile_ccn_number',status='old',
     334     &       form='formatted',iostat=ierr)
     335            if (ierr.eq.0) then
     336              read(91,*) qsurf(iq)
     337              do ilayer=1,nlayermx
     338                read(91,*) q(ilayer,iq)
     339              enddo
     340            else
     341              write(*,*) "No profile_ccn_number file!"
     342            endif
     343            close(91)
     344          endif ! of if (txt.eq."ccn_number")
    314345        enddo ! of do iq=1,nqmx
    315         ! NB: some more initializations (chemistry) is done later
    316346
    317347      else
    … …  
    323353        enddo
    324354      endif ! of if (tracer)
     355     
     356      !write(*,*) "testphys1d q", q(1,:)
     357      !write(*,*) "testphys1d qsurf", qsurf
    325358
    326359c  Date et heure locale du debut du run
    … …  
    613646c       appel de la physique
    614647c       --------------------
     648!      write(*,*) "testphys1d avant q", q(1,:)
    615649      CALL physiq (1,llm,nqmx,
    616650     ,     firstcall,lastcall,
    … …  
    621655C - sorties
    622656     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
     657!      write(*,*) "testphys1d apres q", q(1,:)
     658
    623659
    624660c       evolution du vent : modele 1D

  • trunk/LMDZ.MARS/libf/phymars/tracer.h

    r324 r358  
    1414      real rho_ice     ! Water ice density (kg.m-3)
    1515      real nuice_ref   ! Effective variance of the water ice dist.
     16      real nuice_sed   ! Sedimentation effective variance of the water ice dist.
    1617      real ref_r0        ! for computing reff=ref_r0*r0 (in log.n. distribution)
    1718
    … …  
    2728      integer :: igcm_dust_number ! dust number mixing ratio
    2829                                  !   (transported dust)
     30      integer :: igcm_ccn_mass   ! CCN mass mixing ratio
     31      integer :: igcm_ccn_number ! CCN number mixing ratio
    2932      integer :: igcm_dust_submicron ! submicron dust mixing ratio
    3033                                     !   (transported dust)
    … …  
    6972!     split them in groups (reals, integers and characters)
    7073      COMMON/tracer/radius,rho_q,alpha_lift,alpha_devil,mmol,           &
    71      & varian,r3n_q,rho_dust,rho_ice,nuice_ref,ref_r0,dryness
     74     & varian,r3n_q,rho_dust,rho_ice,nuice_ref,nuice_sed,               &
     75     & ref_r0,dryness
    7276      COMMON/tracer2/                                                   &
    73      & igcm_dustbin,igcm_dust_mass,igcm_dust_number,igcm_dust_submicron,&
     77     & igcm_dustbin,igcm_dust_mass,igcm_dust_number,                    &
     78     & igcm_ccn_mass,igcm_ccn_number,igcm_dust_submicron,               &
    7479     & igcm_h2o_vap,igcm_h2o_ice,igcm_co2,igcm_co,igcm_o,igcm_o1d,      &
    7580     & igcm_o2,igcm_o3,igcm_h,igcm_h2,igcm_oh,igcm_ho2,igcm_h2o2,       &

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

    r38 r358  
    4747      REAL pq(ngrid,nlayer,nqmx)
    4848      real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
    49       real reffdust(ngridmx,nlayermx) ! Dust effective radius (m)
    50       real nueffdust(ngridmx,nlayermx) ! Dust effective variance
    5149
    5250c     Outputs:
    … …  
    7876      EXTERNAL CBRT
    7977
     78      real nueffdust(ngridmx,nlayermx) ! Dust effective variance
     79
    8080c     Local saved variables:
    8181c     ---------------------
    8282
    8383c==================================================================
    84 c 1. Radius used in the physical subroutines
    85 c   (but not in the rad. transfer)
    86 c==================================================================
    8784
    88 c     1.1 Dust particles
    89 c     ------------------
     85      IF (firstcall) THEN
     86c       At firstcall, rdust and rice are not known; therefore
     87c         they need to be computed below.
    9088
    91       IF (doubleq.AND.active) THEN
    92         DO l=1,nlayer
    93           DO ig=1, ngrid
    94             reffdust(ig,l) = ref_r0 *
    95      &        CBRT(r3n_q*pq(ig,l,igcm_dust_mass)/
    96      &        max(pq(ig,l,igcm_dust_number),0.01))
    97             reffdust(ig,l)=min(max(reffdust(ig,l),1.e-10),500.e-6)
    98             nueffdust(ig,l) = exp(varian**2.)-1.
     89c       1.1 Dust particles
     90c       ------------------
     91        IF (doubleq.AND.active) THEN
     92          DO l=1,nlayer
     93            DO ig=1, ngrid
     94              rdust(ig,l) =
     95     &          CBRT(r3n_q*pq(ig,l,igcm_dust_mass)/
     96     &          max(pq(ig,l,igcm_dust_number),0.01))
     97              rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
     98              nueffdust(ig,l) = exp(varian**2.)-1.
     99             ENDDO
     100           ENDDO
     101        ELSE
     102          DO l=1,nlayer
     103            DO ig=1, ngrid
     104              rdust(ig,l) = 0.8E-6
     105              nueffdust(ig,l) = 0.3
     106            ENDDO
    99107          ENDDO
    100         ENDDO
    101       ELSE
    102         DO l=1,nlayer
    103           DO ig=1, ngrid
    104             reffdust(ig,l) = 1.5E-6
    105             nueffdust(ig,l) = 0.3
     108        ENDIF
     109c       1.2 Water-ice particles
     110c       -----------------------
     111        IF (water.AND.activice) THEN
     112          DO l=1,nlayer
     113            DO ig=1,ngrid
     114              rice(ig,l) = max( CBRT(
     115     &          (pq(ig,l,igcm_h2o_ice)/rho_ice +
     116     &          ccn0*(4./3.)*pi*rdust(ig,l)**3.) /
     117     &          (ccn0*4./3.*pi)),rdust(ig,l) )
     118              nuice(ig,l) = nuice_ref
     119            ENDDO
    106120          ENDDO
    107         ENDDO
    108       ENDIF
    109 
    110       DO l=1,nlayer
    111         DO ig=1, ngrid
    112 c         Geometric mean radius = Effective radius / (1+nueff)^5/2
    113           rdust(ig,l) = reffdust(ig,l)/(1.+nueffdust(ig,l))**2.5
    114         ENDDO
    115       ENDDO
    116 
    117 c     1.2 Water-ice particles
    118 c     -----------------------
    119 
    120       IF (firstcall.AND.water.AND.activice) THEN
    121         DO l=1,nlayer
    122           DO ig=1,ngrid
    123             rice(ig,l) = max( CBRT(
    124      &        (pq(ig,l,igcm_h2o_ice)/rho_ice +
    125      &        ccn0*(4./3.)*pi*rdust(ig,l)**3.) /
    126      &        (ccn0*4./3.*pi)),rdust(ig,l) )
    127             nuice(ig,l) = nuice_ref
    128           ENDDO
    129         ENDDO
     121        ENDIF ! of if (water.AND.activice)
    130122        firstcall = .false.
    131       ENDIF
     123      ENDIF ! of if firstcall
    132124
    133125c==================================================================
    … …  
    142134          DO l=1,nlayer
    143135            DO ig=1,ngrid
    144               reffrad(ig,l,iaer) = reffdust(ig,l)
     136              reffrad(ig,l,iaer) = rdust(ig,l) *
     137     &          (1.e0 + nueffdust(ig,l))**2.5
    145138              nueffrad(ig,l,iaer) = nueffdust(ig,l)
    146139            ENDDO
    … …  
    151144          DO l=1,nlayer
    152145            DO ig=1,ngrid
    153               reffrad(ig,l,iaer) = reffdust(ig,l)
     146              reffrad(ig,l,iaer) = rdust(ig,l) * ref_r0
    154147              nueffrad(ig,l,iaer) = nueffdust(ig,l)
    155148            ENDDO
    … …  
    169162          DO l=1,nlayer
    170163            DO ig=1,ngrid
     164c             About reffice, do not confuse the mass mean radius
     165c             (rayon moyen massique) and the number median radius
     166c             (or geometric mean radius, rayon moyen géométrique).
     167c             rice is a mass mean radius, whereas rdust
     168c             is a geometric mean radius:
     169c             number median rad = mass mean rad x exp(-1.5 sigma0^2)
     170c             (Montmessin et al. 2004 paragraph 30). Therefore:
    171171              reffrad(ig,l,iaer)=rice(ig,l)*(1.+nuice_ref)
    172172              nueffrad(ig,l,iaer)=nuice_ref

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

    r334 r358  
    22     &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
    33     &                pq,pdq,pdqcloud,pdqscloud,pdtcloud,
    4      &                nq,naersize,tau,
    5      &                ccn,rdust,rice,nuice,
    6      &                surfdust,surfice)
     4     &                nq,tau,tauscaling,rdust,rice,nuice,
     5     &                rsedcloud,rhocloud)
    76      IMPLICIT NONE
    87
    98c=======================================================================
    10 c     Treatment of saturation of water vapor
     9c  Water-ice cloud formation
     10c 
     11c  Includes two different schemes:
     12c    - A simplified scheme (see simpleclouds.F)
     13c    - An improved microphysical scheme (see improvedclouds.F)
    1114c
     15c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
     16c           J.-B. Madeleine
    1217c
    13 c     Modif de zq si saturation dans l'atmosphere
    14 c     si zq(ig,l)> zqsat(ig,l) ->    zq(ig,l)=zqsat(ig,l)
    15 c     Le test est effectue de bas en haut. L'eau condensee
    16 c    (si saturation) est remise dans la couche en dessous.
    17 c     L'eau condensee dans la couche du bas est deposee a la surface
    18 c       
    19 c    Modification: Franck Montmessin water ice scheme
    20 c                  Francois Forget : change nuclei density & outputs   
    21 c                  Ehouarn Millour: sept.2008, tracers are now handled
    22 c                                   by name (and not fixed index)
    23 c
     18c  2004 - Oct. 2011
    2419c=======================================================================
    2520
    … …  
    3429#include "tracer.h"
    3530#include "comgeomfi.h"
     31#include "dimradmars.h"
    3632
    3733c   Inputs:
    … …  
    3935
    4036      INTEGER ngrid,nlay
     37      integer nq                 ! nombre de traceurs
    4138      REAL ptimestep             ! pas de temps physique (s)
    4239      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    … …  
    5148      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
    5249
    53       integer nq         ! nombre de traceurs
    54       integer naersize   ! nombre de traceurs radiativement actifs (=naerkind)
    55       REAL tau(ngridmx,naersize)
    56       REAL ccn(ngridmx,nlayermx)   ! Cloud condensation nuclei
    57                                    !   (particules kg-1)
     50      REAL tau(ngridmx,naerkind)   ! Column dust optical depth at each point
     51      REAL tauscaling(ngridmx)     ! Convertion factor for dust amount
    5852      real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
    5953
    … …  
    7064      REAL nuice(ngrid,nlay)   ! Estimated effective variance
    7165                               !   of the size distribution
    72 
    73       REAL surfdust(ngrid,nlay)   ! dust surface area (micron2/cm3, if photochemistry)
    74       REAL surfice(ngrid,nlay)    !  ice surface area (micron2/cm3, if photochemistry)
     66      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
     67      real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
    7568
    7669c   local:
    7770c   ------
    7871
    79       REAL CBRT
    80       EXTERNAL CBRT
    8172      INTEGER ig,l
    82 
    83 
    84       REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
    85       REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
    86       REAL zqsat(ngridmx,nlayermx)    ! saturation
    87       REAL zt(ngridmx,nlayermx)       ! local value of temperature
    88 
    89       REAL masse (ngridmx,nlayermx)
    90       REAL epaisseur (ngridmx,nlayermx)
    91       REAL npart(ngridmx,nlayermx)    !  Cloud condensation nuclei (#.cm-3)
    92 !      REAL rfinal        ! Ice crystal radius after condensation(m)
    93       REAL*8 seq           ! Equilibrium saturation ration (accounting for curvature effect)
    94       REAL*8 dzq           ! masse de glace echangee (kg/kg)
    95       REAL lw       !Latent heat of sublimation (J.kg-1)
    96       REAL,PARAMETER :: To=273.15 ! reference temperature, T=273.15 K
    97 
    98       REAL Ctot
    99       REAL*8 ph2o,satu
    100       REAL*8 gr,Cste,up,dwn,newvap
    101 
    10273      LOGICAL,SAVE :: firstcall=.true.
    103 ! To use more refined microphysics, set improved to .true.
    104       LOGICAL,PARAMETER :: improved=.true.
    105 
    106 c     Pour diagnostique :
    107 c     ~~~~~~~~~~~~~~~~~
    108 c     REAL icetot(ngridmx)             ! Total mass of water ice (kg/m2)
    109 c     REAL rave(ngridmx)               ! Mean crystal radius in a column (m)
    110 
    111 ! indexes of water vapour, water ice and dust tracers:
    112       INTEGER,SAVE :: i_h2o=0  ! water vapour
    113       INTEGER,SAVE :: i_ice=0  ! water ice
    11474
    11575c    ** un petit test de coherence
    … …  
    13191        endif
    13292         
    133         i_h2o=igcm_h2o_vap
    134         i_ice=igcm_h2o_ice
    135        
    136         write(*,*) "watercloud: i_h2o=",i_h2o
    137         write(*,*) "            i_ice=",i_ice
     93        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
     94        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
    13895
    13996        firstcall=.false.
    14097      ENDIF ! of IF (firstcall)
    14198
     99c     Main call to the different cloud schemes:
     100      IF (microphys) THEN
     101        CALL improvedclouds(ngrid,nlay,ptimestep,
     102     &             pplay,pt,pdt,
     103     &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     104     &             nq,tauscaling,rdust,rice,nuice,
     105     &             rsedcloud,rhocloud)
     106      ELSE
     107        CALL simpleclouds(ngrid,nlay,ptimestep,
     108     &             pplev,pplay,pzlev,pzlay,pt,pdt,
     109     &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     110     &             nq,tau,rice,nuice,rsedcloud)
     111      ENDIF
    142112
    143 c-----------------------------------------------------------------------
    144 c    1. initialisation
    145 c    -----------------
     113c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
     114c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     115c     Then that should not affect the ice particle radius
     116      do ig=1,ngridmx
     117        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
     118          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
     119     &    rice(ig,2)=rice(ig,3)
     120          rice(ig,1)=rice(ig,2)
     121        end if
     122      end do
    146123
    147 c    On "update" la valeur de q(nqmx) (water vapor) et temperature.
    148 c    On effectue qqes calculs preliminaires sur les couches :
    149 c    masse (kg.m-2), epaisseur(m).
     124c=======================================================================
    150125
    151       do l=1,nlay
    152         do ig=1,ngrid
    153           zq(ig,l,i_h2o)=pq(ig,l,i_h2o)+pdq(ig,l,i_h2o)*ptimestep
    154           zq(ig,l,i_h2o)=max(zq(ig,l,i_h2o),1.E-30) ! FF 12/2004
    155           zq0(ig,l,i_h2o)=zq(ig,l,i_h2o)
    156           zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
    157           masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
    158           epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
     126!!!!!!!!!! FOR PHOTOCHEMISTRY, REIMPLEMENT output of surfdust/surfice
     127!!      if (photochem) then
     128!!c        computation of dust and ice surface area (micron2/cm3)
     129!!c        for heterogeneous chemistry
     130!!
     131!!         do l = 1,nlay
     132!!            do ig = 1,ngrid
     133!!c                                                                     
     134!!c              npart: number density of ccn in #/cm3   
     135!!c                                                     
     136!!               npart(ig,l) = 1.e-6*ccn(ig,l) 
     137!!     $                       *masse(ig,l)/epaisseur(ig,l)
     138!!c
     139!!c              dust and ice surface area
     140!!c                                                                 
     141!!               surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2
     142!!c                                                                 
     143!!               if (rice(ig,l) .ge. rdust(ig,l)) then                   
     144!!                  surfice(ig,l)  = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2
     145!!                  surfdust(ig,l) = 0.
     146!!               else                                                   
     147!!                  surfice(ig,l) = 0.                                 
     148!!               end if                                             
     149!!            end do      ! of do ig=1,ngrid
     150!!         end do         ! of do l=1,nlay
     151!!      end if            ! of photochem
    159152
    160           zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
    161           zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.) ! FF 12/2004
    162           zq0(ig,l,i_ice)=zq(ig,l,i_ice)
    163 
    164 c         This typical profile is not used anymore; rdust is
    165 c           set up in updatereffrad.F.
    166 c         rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
    167         enddo
    168       enddo
    169 
    170       do l=1,nlay
    171         do ig=1,ngrid
    172 c         Calcul du rayon moyen des particules de glace.
    173 c         Hypothese : Dans une couche, la glace presente se
    174 c         repartit uniformement autour du nbre de poussieres
    175 c         dont le rayon moyen est prescrit par rdust.
    176 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    177           rice(ig,l)=CBRT( ( zq(ig,l,i_ice)/rho_ice+
    178      &       ccn(ig,l)*(4./3.)*pi*rdust(ig,l)**3. )
    179      &       / (ccn(ig,l)*4./3.*pi) )
    180           rice(ig,l)=max(rice(ig,l),rdust(ig,l))
    181 c         Effective variance of the size distribution
    182           nuice(ig,l)=nuice_ref
    183         enddo ! of do ig=1,ngrid
    184       enddo ! of dol=1,nlay
    185 
    186       if (photochem) then
    187 c        computation of dust and ice surface area (micron2/cm3)
    188 c        for heterogeneous chemistry
    189 
    190          do l = 1,nlay
    191             do ig = 1,ngrid
    192 c                                                                     
    193 c              npart: number density of ccn in #/cm3   
    194 c                                                     
    195                npart(ig,l) = 1.e-6*ccn(ig,l) 
    196      $                       *masse(ig,l)/epaisseur(ig,l)
    197 c
    198 c              dust and ice surface area
    199 c                                                                 
    200                surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2
    201 c                                                                 
    202                if (rice(ig,l) .ge. rdust(ig,l)) then                   
    203                   surfice(ig,l)  = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2
    204                   surfdust(ig,l) = 0.
    205                else                                                   
    206                   surfice(ig,l) = 0.                                 
    207                end if                                             
    208             end do      ! of do ig=1,ngrid
    209          end do         ! of do l=1,nlay
    210       end if            ! of photochem
    211 c
    212       pdqscloud(1:ngrid,1:nq)=0
    213       pdqcloud(1:ngrid,1:nlay,1:nq)=0
    214       pdtcloud(1:ngrid,1:nlay)=0
    215 
    216 c     icetot(1:ngrid)=0
    217 c     rave(1:ngrid)=0
    218 
    219 c    ----------------------------------------------
    220 c
    221 c
    222 c       Rapport de melange a saturation dans la couche l : -------
    223 c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    224 
    225         call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
    226 
    227 c       taux de condensation (kg/kg/s-1) dans les differentes couches
    228 c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    229 
    230 c       Iceparty is not used anymore: water=>iceparty (JBM).
    231 c       if(iceparty) then
    232 
    233         do l=1,nlay
    234           do ig=1,ngrid
    235 
    236           IF (improved) then
    237 c         Improved microphysics scheme
    238 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    239            Ctot = zq(ig,l,i_h2o) + zq(ig,l,i_ice)
    240            ph2o = zq(ig,l,i_h2o) * 44. / 18. * pplay(ig,l)
    241            satu = zq(ig,l,i_h2o) / zqsat(ig,l)
    242 
    243            call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
    244      &        ph2o,ph2o/satu,seq,rice(ig,l),gr)
    245            Cste = ptimestep * 4. * pi * rice(ig,l)
    246      *              * rho_ice * ccn(ig,l)
    247            up   = zq(ig,l,i_h2o) + Cste * gr * seq
    248            dwn  =    1.       + Cste * gr / zqsat(ig,l)
    249            newvap = min(up/dwn,Ctot)
    250 
    251            gr  = gr * ( newvap/zqsat(ig,l) - seq )
    252            dzq = min( max( Cste * gr,-zq(ig,l,i_ice) )
    253      *             , zq(ig,l,i_h2o) )
    254 
    255 c          Nucleation (sat ratio must be larger than a critical value)
    256 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    257            if (satu.gt.1.) then
    258              if (satu.le.1.4.and.zq(ig,l,i_ice).lt.1.e-8)
    259      *          dzq = 0.
    260            endif
    261 
    262            ELSE
    263 c          Old version
    264 c          ~~~~~~~~~~~
    265            if (zq(ig,l,i_h2o).ge.zqsat(ig,l))then  !  Condensation
    266              dzq=zq(ig,l,i_h2o)-zqsat(ig,l)               
    267            elseif(zq(ig,l,i_h2o).lt.zqsat(ig,l))then  ! Sublimation
    268              dzq=-min(zqsat(ig,l)-zq(ig,l,i_h2o),zq(ig,l,i_ice))
    269            endif
    270 
    271            ENDIF ! of IF (improved)
    272 
    273 c           Water Mass change
    274 c           ~~~~~~~~~~~~~~~~~
    275             zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq
    276             zq(ig,l,i_h2o)=zq(ig,l,i_h2o)-dzq
    277 
    278             rice(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ice
    279      &       +ccn(ig,l)*(4./3.)*pi*rdust(ig,l)**3.)
    280      &       /(ccn(ig,l)*4./3.*pi) ), rdust(ig,l))
    281 
    282             enddo ! of do ig=1,ngrid
    283           enddo ! of do l=1,nlay
    284 
    285 c       The following part have been commented because iceparty
    286 c           is not used anymore: water=>iceparty (JBM).
    287 
    288 c       else   ! if not iceparty
    289 
    290 c         Saturation couche nlay a 2 :
    291 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~
    292 c         do l=nlay,2, -1
    293 c          do ig=1,ngrid
    294 c           if (zq(ig,l,i_h2o).gt.zqsat(ig,l))then
    295 c             zq(ig,l-1,i_h2o)= zq(ig,l-1,i_h2o)+
    296 c    &                          (zq(ig,l,i_h2o)-zqsat(ig,l))
    297 c    &          *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
    298 c             zq(ig,l,i_h2o)=zqsat(ig,l)
    299 c           endif
    300 c          enddo
    301 c         enddo
    302 
    303 c       Saturation couche l=1 si pas iceparty
    304 c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    305 c         do ig=1,ngridmx
    306 c           if (zq(ig,1,i_h2o).gt.zqsat(ig,1))then
    307 c             pdqscloud(ig,i_ice)=(zq(ig,1,i_h2o)-zqsat(ig,1))
    308 c    &           *(pplev(ig,1)-pplev(ig,2))/(g*ptimestep)
    309 c             zq(ig,1,i_h2o)=zqsat(ig,1)
    310 c           endif
    311 c         enddo
    312 
    313 c       endif   ! of if (iceparty)
    314 
    315 c       Tendance finale
    316 c       ~~~~~~~~~~~~~~~
    317         do l=1, nlay
    318           do ig=1,ngridmx
    319             pdqcloud(ig,l,i_h2o)=(zq(ig,l,i_h2o)
    320      &                              -zq0(ig,l,i_h2o))/ptimestep
    321             pdqcloud(ig,l,i_ice) =
    322      &        (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep
    323             lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
    324             pdtcloud(ig,l)=-pdqcloud(ig,l,i_h2o)*lw/cpp
    325           end do
    326         end do
    327 
    328 c       A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
    329 c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    330 c       Then that should not affect the ice particle radius
    331         do ig=1,ngridmx
    332           if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
    333             if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
    334      &      rice(ig,2)=rice(ig,3)
    335             rice(ig,1)=rice(ig,2)
    336           end if
    337         end do
    338        
    339 c**************************************************
    340 c       Output 
    341 c**************************************************
    342 ! NB: for diagnostics use zq(), the updated value of tracers
    343 
    344 c        do ig=1,ngridmx
    345 c         do l=1 ,nlay
    346 c           masse de glace d'eau dans la couche l
    347 c           icetot(ig)=icetot(ig)+masse(ig,l)*zq(ig,l,i_ice)
    348 c           rayon moyen des cristaux dans la colonne ig
    349 c           rave(ig)=rave(ig)+masse(ig,l)*zq(ig,l,i_ice)*rice(ig,l)
    350 c         enddo
    351 c         rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
    352 c         if (icetot(ig)*1000.lt.0.01) rave(ig)=0.
    353 c        enddo   ! (ngridmx)
    354 c**************************************************
    355153
    356154      RETURN
Note: See TracChangeset for help on using the changeset viewer.