Changeset 740 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Jul 25, 2012, 1:03:19 PM (12 years ago)
Author:
tnavarro
Message:

module for ice and radius radius computation

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

Legend:

Unmodified
Added
Removed

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

    r630 r740  
    389389     &                rdust,rice,nuice,
    390390     &                reffrad,nueffrad,
    391      &                pq,tauscaling)
     391     &                pq,tauscaling,tau,pplay)
    392392
    393393c     Computing 3D scattering parameters:

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

    r635 r740  
    33     &                rsedcloud,rhocloud,
    44     &                pq, pdqfi, pdqsed,pdqs_sed,nq,
    5      &                tau, tauscaling)
     5     &                tau,tauscaling)
    66! to use  'getin'
    7       USE ioipsl_getincom
     7      USE ioipsl_getincom
     8      USE updaterad
    89      IMPLICIT NONE
    910
    … …  
    5455c   ------
    5556
    56       REAL CBRT
    57       EXTERNAL CBRT
    58 
    5957      INTEGER l,ig, iq
    6058
    … …  
    6765      real r0dust(ngridmx,nlayermx) ! geometric mean radius used for
    6866                                    !   dust (m)
    69       real r0ccn(ngridmx,nlayermx)  ! geometric mean radius used for
    70                                     !   CCNs (m)
     67!      real r0ccn(ngridmx,nlayermx)  ! geometric mean radius used for
     68!                                    !   CCNs (m)
    7169c     Sedimentation radius of water ice
    7270      real rsedcloud(ngridmx,nlayermx)
    … …  
    186184        ENDIF !of if (doubleq)
    187185
    188         IF (scavenging) THEN
     186        IF (microphys) THEN
    189187          iccn_mass=0
    190188          iccn_number=0
    … …  
    201199            write(*,*) 'callsedim: error! could not identify'
    202200            write(*,*) ' tracers for ccn mass and number mixing'
    203             write(*,*) ' ratio and scavenging is activated!'
     201            write(*,*) ' ratio and microphys is activated!'
    204202            stop
    205203          endif
    206         ENDIF !of if (scavenging)
     204        ENDIF !of if (microphys)
    207205
    208206        IF (water) THEN
    … …  
    255253        do l=1,nlay
    256254          do ig=1, ngrid
    257             r0dust(ig,l) =
    258      &        CBRT(r3n_q*zqi(ig,l,idust_mass)/
    259      &        max(zqi(ig,l,idust_number),0.01))
    260             r0dust(ig,l)=min(max(r0dust(ig,l),1.e-10),500.e-6)
     255     
     256         call updaterdust(zqi(ig,l,igcm_dust_mass),
     257     &                    zqi(ig,l,igcm_dust_number),r0dust(ig,l),
     258     &                    tauscaling(ig))
     259         
    261260          end do
    262261        end do
    263262      endif
    264       if (scavenging) then
    265         do l=1,nlay
    266           do ig=1, ngrid
    267             r0ccn(ig,l) = rsedcloud(ig,l)/(1.+nuice_sed)**4.5
    268           end do
    269         end do
    270       endif
     263
    271264
    272265c =================================================================
    … …  
    279272          if ((doubleq.and.
    280273     &        ((iq.eq.idust_mass).or.
    281      &         (iq.eq.idust_number)))) then !.or.
    282 !     &        (scavenging.and.
    283 !     &        ((iq.eq.iccn_mass).or.
    284 !     &        (iq.eq.iccn_number)))) then
    285 
     274     &         (iq.eq.idust_number)))) then
     275     
    286276c           Computing size distribution:
    287277c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    294284              end do
    295285              sigma0 = varian
    296 c            else             ! ccn
    297 c              do  l=1,nlay
    298 c                do ig=1, ngrid
    299 c                  r0(ig,l)=r0ccn(ig,l)
    300 c                end do
    301 c              end do
    302 c              sigma0 = sqrt(log(1.+nuice_sed))
    303 c            endif
    304286
    305287c        Computing mass mixing ratio for each particle size
    … …  
    358340
    359341          do ir=1,nr
    360 !             IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then
    361               ! Dust sedimentation
     342         
    362343               call newsedim(ngrid,nlay,1,1,ptimestep,
    363344     &         pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
    364345     &         wq,0.5)
    365 !             ELSE
    366 !               ! CCN sedimentation
    367 !               call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep,
    368 !     &         pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir),
    369 !     &         wq,beta)
    370 !             ENDIF
    371346
    372347c            Tendencies
    … …  
    385360c         WATER CYCLE CASE
    386361c -----------------------------------------------------------------
    387 c          else if (water.and.(iq.eq.igcm_h2o_ice)) then
    388362           else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
    389363     &       .or. (iq .eq. igcm_h2o_ice)) then
    … …  
    437411       DO l = 1, nlay
    438412        DO ig=1,ngrid
    439           rdust(ig,l)=
    440      &      CBRT(r3n_q*zqi(ig,l,idust_mass)/
    441      &      max(zqi(ig,l,idust_number),0.01))
    442           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
     413       
     414     
     415         call updaterdust(zqi(ig,l,igcm_dust_mass),
     416     &                    zqi(ig,l,igcm_dust_number),rdust(ig,l),
     417     &                    tauscaling(ig))     
     418
     419         
    443420        ENDDO
    444421       ENDDO
    … …  
    448425c     -------------------------------------
    449426      if (water) then
    450        IF(scavenging) THEN
     427       IF(microphys) THEN
     428       
     429       
    451430        DO l = 1, nlay
    452431          DO ig=1,ngrid
    453             Mo   = zqi(ig,l,igcm_h2o_ice) +
    454      &             zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
    455             No   = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
    456             rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
    457      &                       +zqi(ig,l,iccn_mass)* tauscaling(ig)
    458      &                        / Mo * rho_dust
    459             rhocloud(ig,l) =
    460      &         min(max(rhocloud(ig,l),rho_ice),rho_dust)
    461             rice(ig,l) =
    462      &        ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
    463            if ((Mo.lt.1.e-15) .or. (No.le.1)) rice(ig,l) = 1.e-8
    464 !           print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
    465 !           print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
     432
     433         call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
     434     &    zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
     435     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
    466436           
    467437          ENDDO
    468438        ENDDO
     439       
    469440       ELSE
     441       
    470442        DO l = 1, nlay
    471443          DO ig=1,ngrid
    472             ccntyp =
    473      &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.)
    474             ccntyp = ccntyp /ccn_factor
    475             rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
    476      &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
    477      &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
     444         
     445            call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
     446     &                      tau(ig,1),zlay(ig,l),rice(ig,l))
     447
    478448          ENDDO
    479449        ENDDO
    480        ENDIF ! of IF(scavenging)
     450       ENDIF ! of IF(microphys)
    481451      endif ! of if (water)
    482452

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

    r689 r740  
    55! to use  'getin'
    66      USE ioipsl_getincom
     7      USE updaterad
    78      implicit none
    89     
    … …  
    6970      REAL*8   derf ! Error function
    7071      !external derf
    71 
    72       REAL CBRT
    73       EXTERNAL CBRT
    74 
     72     
    7573      INTEGER ig,l,i
    7674     
    77 
    7875      REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
    7976      REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
    … …  
    104101
    105102      REAL res      ! Resistance growth
    106 
     103     
    107104
    108105c     Parameters of the size discretization
    … …  
    131128     
    132129      LOGICAL test_flag    ! flag for test/debuging outputs
    133       SAVE    test_flag     
     130      SAVE    test_flag   
     131
     132
     133      REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
     134      REAL res_out(ngrid,nlay)
    134135 
    135136
    … …  
    185186
    186187        do i=1,nbin_cld+1
    187     !        rb_cld(i) = log(rb_cld(i)) 
     188!           rb_cld(i) = log(rb_cld(i)) 
    188189            rb_cld(i) = dlog(rb_cld(i))  !! we save that so that it is not computed
    189190                                         !! at each timestep and gridpoint
    … …  
    217218      cste = 4*pi*rho_ice*ptimestep
    218219
     220      res_out(:,:) = 0
     221      rice(:,:) = 1.e-8
     222
    219223c     Initialize the tendencies
    220224      pdqcloud(1:ngrid,1:nlay,1:nq)=0
    … …  
    268272       IF ( satu .ge. 1. ) THEN         ! if there is condensation
    269273
    270 c         Update rdust from last tendencies
    271         rdust(ig,l)=
    272      &     CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
    273      &      (zq(ig,l,igcm_dust_number)+1.e-30))
    274         rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
     274        call updaterccn(zq(ig,l,igcm_dust_mass),
     275     &          zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
    275276
    276277
    … …  
    332333c       Update Dust particles
    333334        zq(ig,l,igcm_dust_mass)   =
    334      &         zq(ig,l,igcm_dust_mass)   - dM/ tauscaling(ig)
     335     &  zq(ig,l,igcm_dust_mass)   - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    335336        zq(ig,l,igcm_dust_number) =
    336      &         zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig)
     337     &  zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    337338c       Update CCNs
    338339        zq(ig,l,igcm_ccn_mass)   =
    339      &         zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig)
     340     &  zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    340341        zq(ig,l,igcm_ccn_number) =
    341      &         zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
     342     &  zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    342343       
    343344        ENDIF ! of is satu >1
    … …  
    353354       IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
    354355
    355 
    356         Mo   = zq(ig,l,igcm_h2o_ice) +
    357      &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
     356     
     357        call updaterice_micro(zq(ig,l,igcm_h2o_ice),
     358     &    zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
     359     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
     360
    358361        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
    359362
    360 
    361         rhocloud(ig,l) =  zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
    362      &                  + zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
    363      &                  / Mo * rho_dust
    364         rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
    365          
    366 
    367 c       ice crystal radius   
    368         rice (ig,l) =
    369      &      CBRT( real(Mo/No) * 0.75 / pi / rhocloud(ig,l) )
    370 
    371 
    372363c       saturation at equilibrium
    373         seq  = exp( 2.*sig(zt(ig,l))*mh2o /
    374      &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )
    375 
     364c       rice should not be too small, otherwise seq value is not valid
     365        seq  = exp(2.*sig(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*
     366     &             max(rice(ig,l),1.e-7)))
     367       
    376368c       get resistance growth
    377369        call growthrate(zt(ig,l),pplay(ig,l),
    378370     &             real(ph2o/satu),rice(ig,l),res)
    379371
     372        res_out(ig,l) = res
    380373
    381374ccccccc  implicit scheme of mass growth
    … …  
    434427     
    435428      ! Get cloud tendencies
    436       pdqcloud(1:ngrid,1:nlay,1:nq) =
    437      & (zq(1:ngrid,1:nlay,1:nq)-zq0(1:ngrid,1:nlay,1:nq))/ptimestep
    438      
    439      
     429        pdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
     430     &   (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
     431     &    zq0(1:ngrid,1:nlay,igcm_h2o_vap))/ptimestep
     432        pdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
     433     &   (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
     434     &    zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
     435        pdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
     436     &   (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
     437     &    zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
     438        pdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
     439     &   (zq(1:ngrid,1:nlay,igcm_ccn_number) -
     440     &    zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
     441     
     442      if (scavenging) then
     443     
     444        pdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
     445     &   (zq(1:ngrid,1:nlay,igcm_dust_mass) -
     446     &    zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
     447        pdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
     448     &   (zq(1:ngrid,1:nlay,igcm_dust_number) -
     449     &    zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
     450         
     451      endif
     452     
     453     
    440454     
    441455!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
    442456!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
    443 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
     457!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
    444458      IF (test_flag) then
    445459     
    446460!       error2d(:) = 0.
    447 !       DO l=1,nlay
    448 !       DO ig=1,ngrid
     461       DO l=1,nlay
     462       DO ig=1,ngrid
    449463!         error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
    450 !       ENDDO
    451 !       ENDDO
     464          satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
     465          satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
     466       ENDDO
     467       ENDDO
    452468
    453469      print*, 'count is ',countcells, ' i.e. ',
    … …  
    470486!         call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
    471487!     &                    error_out)
    472 !         call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1,
    473 !     &                    satubf)
    474 !         call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1,
    475 !     &                    satuaf)
    476 !         call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1,
    477 !     &                    zq0(1,:,igcm_h2o_vap))
    478 !         call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1,
    479 !     &                    zq(1,:,igcm_h2o_vap))
    480 !         call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1,
    481 !     &                    zq0(1,:,igcm_h2o_ice))
    482 !         call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1,
    483 !     &                    zq(1,:,igcm_h2o_ice))
    484 !         call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1,
    485 !     &                    ccnbf)
    486 !         call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1,
    487 !     &                    ccnaf)
     488         call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
     489     &                    res_out)
     490         call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
     491     &                    satubf)
     492         call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
     493     &                    satuaf)
     494         call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
     495     &                    zq0(1,:,igcm_h2o_vap))
     496         call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
     497     &                    zq(1,:,igcm_h2o_vap))
     498         call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
     499     &                    zq0(1,:,igcm_h2o_ice))
     500         call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
     501     &                    zq(1,:,igcm_h2o_ice))
     502         call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
     503     &                    zq0(1,:,igcm_ccn_number))
     504         call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
     505     &                    zq(1,:,igcm_ccn_number))
    488506c         call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
    489507c     &                    gr_out)
    … …  
    494512c         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
    495513c     &                    dN_out)
    496 !         call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
    497 !     &                    zqsat)
     514         call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
     515     &                    zqsat)
    498516!         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
    499517!     &                    satu_out)
    500 !         call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
    501 !     &                    rice)
     518         call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
     519     &                    rice)
    502520!         call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
    503521!     &                    rdust)

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

    r705 r740  
    477477! if scavenging is used, then dustbin should be > 0
    478478
    479          if (scavenging.and.(dustbin.lt.1)) then
    480            print*,'if scavenging is used, then dustbin should > 0'
    481            stop
    482          endif
    483479         if ((microphys.and..not.doubleq).or.
    484      &       (microphys.and..not.active).or.
    485      &       (microphys.and..not.scavenging).or.
    486480     &       (microphys.and..not.water)) then
    487            print*,'if microphys is used, then doubleq,'
    488            print*,'active, water and scavenging must all be used!'
    489            stop
    490          endif
    491          if ((scavenging.and..not.doubleq).or.
    492      &       (scavenging.and..not.active).or.
    493      &       (scavenging.and..not.water).or.
    494      &       (scavenging.and..not.microphys)) then
    495            print*,'if scavenging is used, then doubleq,'
    496            print*,'active, water and microphys must all be used!'
    497            stop
     481             print*,'if microphys is used, then doubleq,'
     482             print*,'and water must be used!'
     483             stop
     484         endif
     485         if (microphys.and..not.scavenging) then
     486             print*,''
     487             print*,'----------------WARNING-----------------'
     488             print*,'microphys is used without scavenging !!!'
     489             print*,'----------------WARNING-----------------'
     490             print*,''
     491         endif
     492
     493         if ((scavenging.and..not.microphys).or.
     494     &       (scavenging.and.(dustbin.lt.1))) then
     495             print*,'if scavenging is used, then microphys'
     496             print*,'must be used!'
     497             stop
    498498         endif
    499499

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

    r648 r740  
    162162        enddo
    163163      endif ! of if (doubleq)
    164       if (scavenging) then
     164      if (microphys) then
    165165        do iq=1,nqmx
    166166          if (noms(iq).eq."ccn_mass") then
    … …  
    173173          endif
    174174        enddo
    175       endif ! of if (scavenging)
     175      endif ! of if (microphys)
    176176      if (submicron) then
    177177        do iq=1,nqmx

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

    r645 r740  
    33     &             pq,pdq,pdqcloud,pdtcloud,
    44     &             nq,tau,rice)
     5      USE updaterad
    56      implicit none
    67c------------------------------------------------------------------
    … …  
    4142      integer nq                 ! nombre de traceurs
    4243      REAL ptimestep             ! pas de temps physique (s)
    43 !      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    4444      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    4545      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
    … …  
    6868      SAVE firstcall
    6969
    70       REAL CBRT
    71       EXTERNAL CBRT
     70           
     71      REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
    7272
    7373      INTEGER ig,l
    … …  
    8686c      REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant
    8787     
    88       REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
    8988
    9089c-----------------------------------------------------------------------
    … …  
    110109      enddo
    111110
    112       do l=1,nlay
    113         do ig=1,ngrid
    114 
    115 c         Typical dust radius profile:
    116           rdusttyp(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
    117 
    118 c         Typical CCN profile:
    119 c         Corrected equation, following Montmessin et al. 2004
    120 c           (N0=2e6 m-3 has been converted into N0=1.3e8 kg-1, otherwise
    121 c            the equation for rice is not homogeneous...)
    122           ccntyp(ig,l)=
    123      &       1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
    124 c         The previously used profile was not correct:
    125 c         ccntyp(ig,l)=( epaisseur(ig,l)/masse(ig,l) ) *
    126 c    &       2.e+6/0.1*max(tau(ig,1),0.001)*exp(-pzlay(ig,l)/10000.)
    127 
    128 c         Reduce number of nuclei
    129 !         TEMPORAIRE : decrease the number of CCNs FF 04/200
    130 !         reduction facteur 3
    131 !         ccntyp(ig,l) = ccntyp(ig,l) / 27.
    132 !         reduction facteur 2
    133 !         ccntyp(ig,l) = ccntyp(ig,l) / 8.
    134 c -----------------------------------------------------------------
    135           ccntyp(ig,l) = ccntyp(ig,l) / ccn_factor
    136 
    137         enddo ! of do ig=1,ngrid
    138       enddo ! of do l=1,nlay
    139111
    140112      pdqcloud(1:ngrid,1:nlay,1:nq)=0
    … …  
    167139          zq(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)-dzq
    168140         
    169           Mcon_out(ig,l) = dzq
    170 
    171 c         Calcul du rayon moyen des particules de glace.
    172 c         Hypothese : Dans une couche, la glace presente se
    173 c         repartit uniformement autour du nbre de poussieres
    174 c         dont le rayon moyen est prescrit par rdusttyp.
    175 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    176           rice(ig,l)=max( CBRT ( (zq(ig,l,igcm_h2o_ice)/rho_ice
    177      &      +ccntyp(ig,l)*(4./3.)*pi*rdusttyp(ig,l)**3.)
    178      &      /(ccntyp(ig,l)*4./3.*pi) ), rdusttyp(ig,l))
    179 
    180141
    181142        enddo ! of do ig=1,ngrid
    … …  
    195156      end do
    196157
    197 c------------------------------------------------------------------
    198 c     TEST_JBM
    199 !      IF (ngrid.eq.1) THEN
    200 !         call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
    201 !     &                    Mcon_out)
    202 !         call WRITEDIAGFI(ngrid,"rdusttyp","rdusttyp","m",1,
    203 !     &                    rdusttyp)
    204 !         call WRITEDIAGFI(ngrid,"ccntyp","ccntyp","kg-1",1,
    205 !     &                    ccntyp)
    206 !      ENDIF
     158c     ice crystal radius
     159      do l=1, nlay
     160        do ig=1,ngridmx
     161          call updaterice_typ(zq(ig,l,igcm_h2o_ice),
     162     &       tau(ig,1),pzlay(ig,l),rice(ig,l))
     163        end do
     164      end do
     165
    207166c------------------------------------------------------------------
    208167      return

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

    r712 r740  
    307307          longwatercaptag(203)    = .true. ! outlier, lat = 75
    308308          longwatercaptag(207)    = .true. ! outlier, lat = 75
    309           longwatercaptag(242)    = .true. ! outlier, lat = 75
    310309          longwatercaptag(244)    = .true. ! outlier, lat = 75
    311310          longwatercaptag(246)    = .true. ! outlier, lat = 75
    312           longwatercaptag(248)    = .true. ! outlier, lat = 75
    313311          longwatercaptag(250)    = .true. ! outlier, lat = 75
    314312          longwatercaptag(252)    = .true. ! outlier, lat = 75
    315313          longwatercaptag(254)    = .true. ! outlier, lat = 75
    316           longwatercaptag(256:257)= .true. ! outlier, lat = 75
     314          longwatercaptag(256)    = .true. ! outlier, lat = 75
    317315          endif
    318316!--------------------------------------------------------------       

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

    r703 r740  
    66! to use  'getin'
    77      USE ioipsl_getincom
     8      USE updaterad
    89      IMPLICIT NONE
    910
    … …  
    9293      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
    9394      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
    94      
    95       REAL CBRT
    96       EXTERNAL CBRT
    9795
    9896
    … …  
    261259          ENDDO
    262260       ENDDO
     261       
    263262c------ Tracers tendencies pdqcloud
    264263       DO l=1,nlay
    … …  
    270269     &        subpdq(ig,l,igcm_h2o_vap)/real(imicro)
    271270     &       - pdq(ig,l,igcm_h2o_vap)
    272           IF(microphys) THEN
    273             pdqcloud(ig,l,igcm_dust_mass) =
    274      &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
    275      &       - pdq(ig,l,igcm_dust_mass)
    276             pdqcloud(ig,l,igcm_dust_number) =
    277      &        subpdq(ig,l,igcm_dust_number)/real(imicro)
    278      &       - pdq(ig,l,igcm_dust_number)
     271         ENDDO
     272       ENDDO
     273       
     274       IF(microphys) THEN
     275        DO l=1,nlay
     276         DO ig=1,ngrid
    279277            pdqcloud(ig,l,igcm_ccn_mass) =
    280278     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro)
    … …  
    283281     &        subpdq(ig,l,igcm_ccn_number)/real(imicro)
    284282     &       - pdq(ig,l,igcm_ccn_number)
    285            ENDIF
    286          ENDDO
    287        ENDDO
     283         ENDDO
     284        ENDDO
     285       ENDIF
     286       
     287       IF(scavenging) THEN
     288        DO l=1,nlay
     289         DO ig=1,ngrid
     290            pdqcloud(ig,l,igcm_dust_mass) =
     291     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
     292     &       - pdq(ig,l,igcm_dust_mass)
     293            pdqcloud(ig,l,igcm_dust_number) =
     294     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
     295     &       - pdq(ig,l,igcm_dust_number)
     296         ENDDO
     297        ENDDO
     298       ENDIF
    288299
    289300c------- Due to stepped entry, other processes tendencies can add up to negative values
    290301c------- Therefore, enforce positive values and conserve mass
    291302
     303
    292304       IF(microphys) THEN
     305        DO l=1,nlay
     306         DO ig=1,ngrid
     307          IF ((pq(ig,l,igcm_ccn_number) +
     308     &      ptimestep* (pdq(ig,l,igcm_ccn_number) +
     309     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
     310     &   .or. (pq(ig,l,igcm_ccn_mass) +
     311     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) +
     312     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
     313         pdqcloud(ig,l,igcm_ccn_number) =
     314     &     - pq(ig,l,igcm_ccn_number)/ptimestep
     315     &     - pdq(ig,l,igcm_ccn_number) + 1.
     316         pdqcloud(ig,l,igcm_dust_number) = 
     317     &     -pdqcloud(ig,l,igcm_ccn_number)
     318         pdqcloud(ig,l,igcm_ccn_mass) =
     319     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
     320     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
     321         pdqcloud(ig,l,igcm_dust_mass) =
     322     &     -pdqcloud(ig,l,igcm_ccn_mass)
     323          ENDIF
     324         ENDDO
     325        ENDDO
     326       ENDIF
     327
     328       IF(scavenging) THEN
    293329        DO l=1,nlay
    294330         DO ig=1,ngrid
    … …  
    310346     &     -pdqcloud(ig,l,igcm_dust_mass)
    311347          ENDIF
    312           IF ((pq(ig,l,igcm_ccn_number) +
    313      &      ptimestep* (pdq(ig,l,igcm_ccn_number) +
    314      &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
    315      &   .or. (pq(ig,l,igcm_ccn_mass) +
    316      &      ptimestep* (pdq(ig,l,igcm_ccn_mass) +
    317      &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
    318          pdqcloud(ig,l,igcm_ccn_number) =
    319      &     - pq(ig,l,igcm_ccn_number)/ptimestep
    320      &     - pdq(ig,l,igcm_ccn_number) + 1.
    321          pdqcloud(ig,l,igcm_dust_number) = 
    322      &     -pdqcloud(ig,l,igcm_ccn_number)
    323          pdqcloud(ig,l,igcm_ccn_mass) =
    324      &     - pq(ig,l,igcm_ccn_mass)/ptimestep
    325      &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
    326          pdqcloud(ig,l,igcm_dust_mass) =
    327      &     -pdqcloud(ig,l,igcm_ccn_mass)
    328           ENDIF
    329348         ENDDO
    330349        ENDDO
    … …  
    358377         DO ig=1,ngrid
    359378
    360           rdust(ig,l)=
    361      &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
    362      &           (pdq(ig,l,igcm_dust_mass)
    363      &           +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
    364      &      (pq(ig,l,igcm_dust_number)+
    365      &           (pdq(ig,l,igcm_dust_number)+
    366      &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
    367           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
    368 
    369           Mo = pq(ig,l,igcm_h2o_ice)
    370      &         + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
    371      &         + (pq(ig,l,igcm_ccn_mass)
    372      &          + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
    373      &         *tauscaling(ig) + 1.e-30
    374           No = (pq(ig,l,igcm_ccn_number)
    375      &          + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
    376      &         *tauscaling(ig) + 1.e-30
    377           rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) + 
    378      &                   pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
    379      &                   / Mo * rho_ice
    380      &                  + (pq(ig,l,igcm_ccn_mass)
    381      &                  + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
    382      &                   *tauscaling(ig)/ Mo * rho_dust
    383           rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
    384           if ((Mo.lt.1.e-15) .or. (No.le.1.)) then
    385               rice(ig,l) = 1.e-8
    386           else
    387               !! AS: only perform computations if conditions not met
    388               rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
    389               rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
    390           endif
    391           rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
    392      &                         (1.+nuice_sed)*(1.+nuice_sed),
    393      &                         rdust(ig,l) )
    394           rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
    395          ENDDO
    396         ENDDO
    397        ELSE
    398        
    399         if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0)) then
    400          ! recompute rdust(), if we have dust_mass & dust_number tracers
    401          DO l=1,nlay
     379        call updaterdust(
     380     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
     381     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
     382     &    pdqcloud(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
     383     &    pq(ig,l,igcm_dust_number) +                 ! dust number
     384     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
     385     &    pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
     386     &    rdust(ig,l))
     387
     388         ENDDO
     389        ENDDO
     390      ENDIF
     391       
     392       
     393      IF(microphys) THEN
     394       
     395       DO l=1, nlay
     396         DO ig=1,ngrid
     397
     398        call updaterice_micro(
     399     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
     400     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
     401     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
     402     &    pq(ig,l,igcm_ccn_mass) +                   ! ccn mass
     403     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
     404     &    pdqcloud(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
     405     &    pq(ig,l,igcm_ccn_number) +                 ! ccn number
     406     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
     407     &    pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
     408     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
     409         
     410         ENDDO
     411       ENDDO
     412       
     413      ELSE ! no microphys
     414       
     415        DO l=1,nlay
    402416          DO ig=1,ngrid
    403            rdust(ig,l)=
    404      &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
    405      &           (pdq(ig,l,igcm_dust_mass)
    406      &           +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
    407      &      (pq(ig,l,igcm_dust_number)+
    408      &           (pdq(ig,l,igcm_dust_number)+
    409      &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
    410            rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
    411           ENDDO
    412          ENDDO
    413         endif ! of if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0))
    414        
    415         DO l=1,nlay
    416           DO ig=1,ngrid
    417             ccntyp =
    418      &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
    419             ccntyp = ccntyp /ccn_factor
    420             rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
    421      &       + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
    422      &      +ccntyp*(4./3.)*pi*rdust(ig,l)*rdust(ig,l)*rdust(ig,l))
    423      &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
    424             rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
    425      &                         (1.+nuice_sed)*(1.+nuice_sed),
    426      &                         rdust(ig,l) )
    427             rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
    428           ENDDO
    429         ENDDO
    430        ENDIF ! of IF(scavenging)
    431        
    432        
    433 ! used for rad. transfer calculations
    434 ! nuice is constant because a lognormal distribution is prescribed
    435       nuice(1:ngrid,1:nlay)=nuice_ref
    436 
    437 
    438 !--------------------------------------------------------------
    439 !--------------------------------------------------------------
    440      
    441 
     417         
     418        call updaterice_typ(
     419     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
     420     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
     421     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
     422     &    tau(ig,l),pzlay(ig,l),rice(ig,l))
     423
     424          ENDDO
     425         ENDDO
     426       
     427       ENDIF ! of IF(microphys)
     428     
     429     
     430     
    442431c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
    443432c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    450439        end if
    451440      end do
     441       
     442       
     443       DO l=1,nlay
     444         DO ig=1,ngrid
     445           rsedcloud(ig,l)=max(rice(ig,l)*
     446     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
     447     &                    rdust(ig,l))
     448!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
     449         ENDDO
     450       ENDDO
     451       
     452! used for rad. transfer calculations
     453! nuice is constant because a lognormal distribution is prescribed
     454      nuice(1:ngrid,1:nlay)=nuice_ref
     455
     456
    452457
    453458c=======================================================================
Note: See TracChangeset for help on using the changeset viewer.