Ignore:
Timestamp:
Jan 22, 2019, 5:02:23 PM (6 years ago)
Author:
mvals
Message:

Mars GCM:

  • Update of rocketduststorm_mod.F90 :

We want to separate both parametrizations related to the formation of the detached dust layers. Therefore, rocketduststorm_mod.F90 now only comprises the rocket dust storm scheme, whereas it contained
also before the calculation of the vertical velocity induced by the presence of the sub-grid scale topography. This latter part is under development and will be integrated as a fully independant
parametrization: the aim is to simulate the entrainment of dust by slope winds, from the boundary layer up to the top of the sub-grid scale topography.

  • Addition of initial surface parameters "summit" and "base" to prepare the previously described slope wind parametrization

MV

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90

    r1985 r2079  
    44
    55      REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1)
    6       REAL, SAVE, ALLOCATABLE :: alpha_hmons(:) ! slope winds lifting mesh fraction
    76     
    87      CONTAINS
    … …  
    1110! ROCKET DUST STORM - vertical transport and detrainment
    1211!=======================================================================
    13 ! calculating the vertical flux
    14 ! calling vl_storm : transport scheme of stormdust tracers
    15 ! detrainement of stormdust into nomal background dust
     12! calculation of the vertical flux
     13! call of vl_storm : Van Leer transport scheme of the dust tracers
     14! detrainement of stormdust in the background dust
    1615! -----------------------------------------------------------------------
    17 ! Authors: C. Wang; F. Forget; T. Bertrand
     16! Authors: M. Vals; C. Wang; F. Forget; T. Bertrand
    1817! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France
    1918! -----------------------------------------------------------------------
    … …  
    2726!             input sub-grid scale cloud
    2827                                 clearsky,totcloudfrac,                &
    29 !                   output
    30                                  pdqrds,wspeed,dsodust,dsords)
    31 
    32       use tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number &
     28!             output
     29                                 pdqrds,wrad,dsodust,dsords)
     30
     31      USE tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number &
    3332                            ,igcm_dust_mass,igcm_dust_number          &
    3433                            ,rho_dust
    3534      USE comcstfi_h, only: r,g,cpp,rcp
    36       use dimradmars_mod, only: albedo,naerkind
    37       use comsaison_h, only: dist_sol,mu0,fract
    38       use surfdat_h, only: emis,co2ice,zmea, zstd, zsig, hmons
    39       use planetwide_mod, only: planetwide_maxval,planetwide_minval
    40 !      use rocketstorm_h, only: rdsinject
    41       use callradite_mod
    42       implicit none
     35      USE dimradmars_mod, only: albedo,naerkind
     36      USE comsaison_h, only: dist_sol,mu0,fract
     37      USE surfdat_h, only: emis,co2ice,zmea, zstd, zsig, hmons
     38      USE callradite_mod
     39      IMPLICIT NONE
    4340
    4441!--------------------------------------------------------
    … …  
    6663
    6764!     input for second radiative transfer
    68       logical, INTENT(IN) :: clearatm
     65      LOGICAL, INTENT(IN) :: clearatm
    6966      INTEGER, INTENT(INOUT) :: icount
    70       real, intent(in) :: zday
    71       real, intent(in) :: zls
    72       real, intent(in) :: tsurf(ngrid)
    73       integer, intent(in) :: igout
    74       real, intent(in) :: totstormfract(ngrid)
     67      REAL, INTENT(IN) :: zday
     68      REAL, INTENT(IN) :: zls
     69      REAL, INTENT(IN) :: tsurf(ngrid)
     70      INTEGER, INTENT(IN) :: igout
     71      REAL, INTENT(IN) :: totstormfract(ngrid)
     72     
    7573!     sbgrid scale water ice clouds
    7674      logical, intent(in) :: clearsky
    … …  
    8280
    8381      REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining
    84       REAL, INTENT(OUT) :: wspeed(ngrid,nlayer+1)   ! vertical speed within the rocket dust storm
     82      REAL, INTENT(OUT) :: wrad(ngrid,nlayer+1)   ! vertical speed within the rocket dust storm
    8583      REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust
    8684      REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust
    … …  
    9088!--------------------------------------------------------
    9189      INTEGER l,ig,tsub,iq,ll
    92 ! chao local variables from callradite.F       
     90!     local variables from callradite.F       
    9391      REAL zdtlw1(ngrid,nlayer)    ! (K/s) storm
    9492      REAL zdtsw1(ngrid,nlayer)    ! (K/s) storm
    … …  
    9795      REAL ztlev(nlayer)     ! temperature at intermediate levels l+1/2 without last level
    9896
    99       REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2
    100       REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer)   ! rad. heating rate at intermediate levels l+1/2
    101 
    102       REAL zqstorm_mass(ngrid,nlayer) ! tracer pq mass intermediate
    103       REAL zqstorm_mass_col(nlayer)
    104       REAL zqstorm_number(ngrid,nlayer) ! tracer field pq number intermediate
    105       REAL zqstorm_number_col(nlayer)
    106 
    107       REAL zqi_mass(ngrid,nlayer)           ! tracer pq mass intermediate
    108       REAL zqi_number(ngrid,nlayer)         ! tracer pq mass intermediate
    109       REAL zdqvlstorm_mass(ngrid,nlayer)    ! tendancy pdq mass after vertical transport
    110       REAL zdqvlstorm_number(ngrid,nlayer)  ! tendancy pdq number after vertical transport
    111       REAL zdqdetstorm_mass(ngrid,nlayer)   ! tendancy field pq mass after detrainment only
    112       REAL zdqdetstorm_number(ngrid,nlayer) ! tendancy field pq number after detrainment only
    113 
    114       REAL zdqenv_mass(ngrid,nlayer)    ! tendancy pdq mass from dust->
    115                                         !  stormdust in slp
    116       REAL zdqenv_number(ngrid,nlayer)  ! tendancy pdq number from dust->
    117                                         !  stormdust in slp
    118 
    119       REAL masse(nlayer)
     97      REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 for stormdust
     98      REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer)   ! rad. heating rate at intermediate levels l+1/2 for background dust
     99
     100      REAL zq_stormdust_mass(ngrid,nlayer) ! intermediate tracer stormdust mass
     101      REAL zq_stormdust_number(ngrid,nlayer) ! intermediate tracer stormdust number
     102      REAL zq_dust_mass(ngrid,nlayer)           ! intermediate tracer dust mass
     103      REAL zq_dust_number(ngrid,nlayer)         ! intermediate tracer dust number
     104
     105      REAL mr_stormdust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust mass)
     106      REAL mr_stormdust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust number)
     107      REAL mr_dust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (dust mass)
     108      REAL mr_dust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (sdust number)
     109
     110      REAL zq_vl_col(nlayer) ! column intermediate tracer used by Van Leer (number)
     111      REAL zn_vl_col(nlayer) ! column intermediate tracer used by Van Leer (mass)
     112                   
     113      REAL dqvl_stormdust_mass(ngrid,nlayer)    ! tendancy of vertical transport (stormdust mass)
     114      REAL dqvl_stormdust_number(ngrid,nlayer)  ! tendancy of vertical transport (stormdust number)
     115      REAL dqvl_dust_mass(ngrid,nlayer)    ! tendancy of vertical transport (dust mass)
     116      REAL dqvl_dust_number(ngrid,nlayer)  ! tendancy of vertical transport (dust number)
     117      REAL dqdet_stormdust_mass(ngrid,nlayer)   ! tendancy of detrainement (stormdust mass)
     118      REAL dqdet_stormdust_number(ngrid,nlayer) ! tendancy of detrainement (stormdust number)
     119
     120      REAL masse(nlayer)     ! mass of atmosphere (kg/m2)
    120121      REAL zq(ngrid,nlayer,nq)   ! updated tracers
    121122     
    122       REAL wrds(nlayer)          ! vertical flux within the rocket dust storm
    123       REAL wqrdsmass(nlayer+1)   ! flux mass from vl_storm
    124       REAL wqrdsnumber(nlayer+1) ! flux number from vl_storm
    125 
    126       INTEGER nsubtimestep    !number of subtimestep when calling vl_storm
    127       REAL subtimestep        !ptimestep/nsubtimestep
    128       REAL dtmax              !considered time needed for dust to cross one layer
    129                               !minimal value over a column
    130       logical storm(ngrid)    !logical : true if you have some storm dust in the column
    131 !     real slope(ngrid)    !logical : true if you don't have storm and have
    132                               !a slope
    133 !     real wslplev(ngrid,nlayer)
    134 !     real wslp(ngrid,nlayer)
    135       REAL coefdetrain           !coefficient for detrainment : % of stormdust detrained
    136 
    137       REAL,PARAMETER:: coefmin =0.025 !C 0<c<1 Minimum fraction of stormdust detrained 
    138       REAL,PARAMETER:: detrainlim =0.1!0.25 !L  stormdust detrained if wspeed < detrainlim 
    139       REAL,PARAMETER:: wlim =10.   ! maximum vertical speed of rocket storms (m/s)
    140       REAL,PARAMETER:: secu=3.      !coefficient on wspeed to avoid dust crossing many layers during subtimestep
    141  
    142       ! terms for buoyancy and W^2 in equation:
    143       !   w*dw/dz = k1*g*(T'-T)/T - k2*w^2
    144       real,parameter:: k1=0.00001
    145       real,parameter:: k2=0.25
    146       real,parameter:: mu0lim=0.1
    147      
    148       ! diagnose
    149       REAL detrainment(ngrid,nlayer)
    150       real lapserate(ngrid,nlayer)
    151       real deltahr(ngrid,nlayer+1)
    152       real stormdust_m0(ngrid,nlayer)
    153       real stormdust_m1(ngrid,nlayer)
    154       real stormdust_m2(ngrid,nlayer)
    155    
    156       real hmax,hmin
    157 !     real zh(ngrid,nlayer)
     123      REAL w(nlayer)          ! air mass flux (calculated with the vertical wind velocity profile) used as input in Van Leer (kgair/m2)
     124      REAL wqmass(nlayer+1)   ! tracer (dust_mass) mass flux in Van Leer (kg/m2)
     125      REAL wqnumber(nlayer+1) ! tracer (dust_number) mass flux in Van Leer (kg/m2)
     126
     127      LOGICAL storm(ngrid)    ! true when there is a dust storm (if the opacity is high): trigger the rocket dust storm scheme
     128      REAL coefdetrain        ! coefficient for detrainment : % of stormdust detrained
     129      INTEGER scheme(ngrid)   ! triggered scheme
     130           
     131      REAL,PARAMETER:: coefmin =0.025 ! 0<coefmin<1 Minimum fraction of stormdust detrained 
     132      REAL,PARAMETER:: wmin =0.1!0.25 ! stormdust detrainment if wrad < wmin 
     133      REAL,PARAMETER:: wmax =10.   ! maximum vertical velocity of the rocket dust storms (m/s)
     134     
     135!     diagnostics
     136      REAL lapserate(ngrid,nlayer)
     137      REAL deltahr(ngrid,nlayer+1)
    158138   
    159       logical,save :: firstcall=.true.
    160       real alpha(ngrid) ! scale of the vertical motion (applicable for
    161                         ! rds and slp
    162       ! variables for radiative transfer
     139      LOGICAL,SAVE :: firstcall=.true.
     140
     141!     variables for the radiative transfer
    163142      REAL  fluxsurf_lw1(ngrid)
    164143      REAL  fluxsurf_sw1(ngrid,2)
    … …  
    175154      REAL  nuice(ngrid,nlayer)
    176155
    177       !variables related to slope,reference layer...
    178 !     integer lref(ngrid),llref ! the reference layer of slopewind
    179 !     real buoyt(nlayer) ! buoyancy term when there is a subgrid mountain
    180       real slpbg(ngrid)  ! temperature difference at half height of a mountain
    181 
    182       real zqdustslp(ngrid,nlayer),zndustslp(ngrid,nlayer)
    183       real zqstormdustslp(ngrid,nlayer),znstormdustslp(ngrid,nlayer)
    184 
    185       real rdsdustqvl0(ngrid,nlayer)
    186       real rdsdustqvl1(ngrid,nlayer)
    187 !     real q2rds(ngrid,nlayer)
    188 
    189       !for second formule of wslp
    190       real wtemp(ngrid,nlayer) ! a intermediate variable for wspeed
    191  
    192       !merge rds and slp
    193       real newzt(ngrid,nlayer) !temperature with perturbation (integrated from
    194                                !  vetical motion)
    195       real w0(ngrid) !prescribed slope winds at the first layer of atmosphere
    196 !     real w1(ngrid) !prescribed slope winds at the first layer of atmosphere
    197 !     real ztb1(ngrid)
    198       real wadiabatic(ngrid,nlayer) !for diagnosis
    199 !     real czt(nlayer),czlay(nlayer),czlev(nlayer+1) !temporary variables
    200       real tnew(nlayer) !interpolated temperature profile above the top of Mons
    201       real envtemp(nlayer) !interpolated background temperature profile
    202                            ! as the same height as tnew
    203       real envt(ngrid,nlayer) ! output,for diagnosing
    204       integer scheme(ngrid)   ! diagnose
    205 
    206       ! Chao: for checking conservation of dust
    207 !     real totdust0(ngrid)
    208 !     real totdust1(ngrid)
    209 
    210 
    211       ! **********************************************************************
    212       ! **********************************************************************
    213       ! Detached dust layers parametrization: two processes are included
    214       ! 1) rocket dust storm
     156
     157      ! **********************************************************************
     158      ! **********************************************************************
     159      ! Rocket dust storm parametrization to reproduce the detached dust layers
     160      ! during the dust storm season:
    215161      !     The radiative warming due to the presence of storm dust is
    216162      !     balanced by the adiabatic cooling. The tracer "storm dust" 
    217163      !     is transported by the upward/downward flow.
    218       ! 2) daytime slope winds
    219       !     The daytime thermally driven upslope wind blows dust from the
    220       !     bottom to the top of the mountain, upward flow keeps rising untill
    221       !     the velocity becomes zero. Both the storm dust and environment dust
    222       !     will be transported.
    223164      ! **********************************************************************
    224165      ! **********************************************************************     
    225166      !!    1. Radiative transfer in storm dust
    226167      !!    2. Compute vertical velocity for storm dust
    227       !!      case 1 storm = false and nighttime: nothing to do
    228       !!      case 2 daytime slope wind scheme: (mu0(ig) > mu0lim and with storm=false)
    229       !!      case 3 rocket dust storm (storm=true)
    230       !!    3. Vertical transport
     168      !!      case 1 storm = false: nothing to do
     169      !!      case 2 rocket dust storm (storm=true)
     170      !!    3. Vertical transport (Van Leer)
    231171      !!    4. Detrainment
    232172      ! **********************************************************************
    233173      ! **********************************************************************
    234      
    235 
    236       !! **********************************************************************
    237       !! Firstcall: Evaluate slope wind mesh fraction 
    238       IF (firstcall) then   
    239         call planetwide_maxval(hmons,hmax )
    240         call planetwide_minval(hmons,hmin )
    241         do ig=1,ngrid
    242           ! It's hard to know the exact the scale of upward flow,
    243           !  we assume that the maximun is 10% of the mesh area.
    244           alpha_hmons(ig)= 0.1*(hmons(ig)-hmin)/(hmax-hmin)
    245         enddo
    246         firstcall = .false.
    247       ENDIF !firstcall
    248      
    249       ! **********************************************************************
    250       ! 0. Initializations
     174
     175     
     176      ! **********************************************************************
     177      ! Initializations
    251178      ! **********************************************************************
    252179      storm(:)=.false.
    253180      pdqrds(:,:,:) = 0.
    254       zdqdetstorm_mass(:,:)=0.
    255       zdqdetstorm_number(:,:)=0.
    256       wspeed(:,:)=0.
    257       detrainment(:,:)=0.
    258       zqstorm_mass_col(:)=0.
    259       zqstorm_number_col(:)=0.
     181      mr_dust_mass(:,:)=0.
     182      mr_dust_number(:,:)=0.
     183      mr_stormdust_mass(:,:)=0.
     184      mr_stormdust_number(:,:)=0.
     185      dqvl_dust_mass(:,:)=0.
     186      dqvl_dust_number(:,:)=0.
     187      dqvl_stormdust_mass(:,:)=0.
     188      dqvl_stormdust_number(:,:)=0.
     189      dqdet_stormdust_mass(:,:)=0.
     190      dqdet_stormdust_number(:,:)=0.
     191      wrad(:,:)=0.
    260192      lapserate(:,:)=0.
    261193      deltahr(:,:)=0.
    262       rdsdustqvl0(:,:)=0.
    263       rdsdustqvl1(:,:)=0.
    264       zqstormdustslp(:,:)=0.
    265       znstormdustslp(:,:)=0.
    266       zqdustslp(:,:)=0.
    267       zndustslp(:,:)=0.
    268       zq(:,:,:) = 0.
    269      
    270       w0(:)=0.
    271 !     w1(:)=0.
    272 !     ztb1(:)=0.
    273       newzt(:,:)=0
    274       wtemp(:,:)=0.
    275       wadiabatic(:,:)=0.
    276       slpbg(:)=0.
    277 !     buoyt(:)=0.
    278       tnew(:)=0.
    279       envtemp(:)=0.
    280       envt(:,:)=0.
    281194      scheme(:)=0
    282       alpha(:)=0.
    283       stormdust_m0(:,:)=0.
    284       stormdust_m1(:,:)=0.
    285       stormdust_m2(:,:)=0.
    286 !     totdust0(:)=0.
    287 !     totdust1(:)=0.
    288195     
    289196      !! no update for the stormdust tracer and temperature fields
    … …  
    292199      zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)
    293200
    294       !! get actual q for stormdust and dust tracers
    295       do l=1,nlayer
    296            do ig=1, ngrid     
    297              zqi_mass(ig,l)=zq(ig,l,igcm_dust_mass)
    298              zqi_number(ig,l)=zq(ig,l,igcm_dust_number)
    299              zqstorm_mass(ig,l)=zq(ig,l,igcm_stormdust_mass)
    300              zqstorm_number(ig,l)=zq(ig,l,igcm_stormdust_number)     
    301              !for diagnostics:
    302              stormdust_m0(ig,l)=zqstorm_mass(ig,l)
    303            enddo
    304       enddo ! of do l=1,nlayer
    305 
    306       !! Check if there is a rocket dust storm
    307       do ig=1,ngrid
     201
     202      zq_dust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_mass)
     203      zq_dust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_number)
     204      zq_stormdust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_mass)
     205      zq_stormdust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_number)
     206
     207      ! *********************************************************************
     208      ! 0. Check if there is a storm
     209      ! *********************************************************************
     210      DO ig=1,ngrid
    308211        storm(ig)=.false.
    309         do l=1,nlayer
    310           if (zqstorm_mass(ig,l)/zqi_mass(ig,l) .gt. 1.E-4) then
     212        DO l=1,nlayer
     213          IF (zq(ig,l,igcm_stormdust_mass) &
     214          .gt. zq(ig,l,igcm_dust_mass)*(1.E-4)) THEN
    311215            storm(ig)=.true.
    312             exit
    313           endif
    314         enddo       
    315       enddo
     216            EXIT
     217          ENDIF
     218        ENDDO     
     219      ENDDO
    316220     
    317221      ! *********************************************************************
    … …  
    328232      ! **********************************************************************
    329233      ! 2. Compute vertical velocity for storm dust
    330       ! **********************************************************************     
    331       DO ig=1,ngrid
    332         !! **********************************************************************
    333         !! 2.1 case 1: Nothing to do when no storm and no slope, or
    334         !!             no storm and not daytime
    335         IF (((alpha_hmons(ig) .EQ. 0.) .AND. .NOT.(storm(ig))) &
    336              .OR. ((mu0(ig) .LE. mu0lim) .AND. .NOT.(storm(ig))) ) then
    337           scheme(ig)=1
    338           cycle
    339         endif
     234      ! **********************************************************************
     235        !! **********************************************************************
     236        !! 2.1 Nothing to do when no storm
     237        !!             no storm
     238        DO ig=1,ngrid       
     239          IF (.NOT.(storm(ig))) then
     240            scheme(ig)=1
     241            cycle
     242          ENDIF ! IF (storm(ig))
     243        ENDDO ! DO ig=1,ngrid                 
    340244       
     245        !! **********************************************************************
     246        !! 2.2 Calculation of the extra heating : computing heating rates
     247        !! gradient at boundaries of each layer, start from surface
     248        DO ig=1,ngrid
     249          zdtvert(1)=0. !This is the env. lapse rate
     250          DO l=1,nlayer-1
     251            zdtvert(l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
     252            lapserate(ig,l+1)=zdtvert(l+1) !for diagnosing
     253          ENDDO
     254     
     255          ! computing heating rates gradient at boundraies of each layer
     256          ! start from surface
     257          zdtlw1_lev(1)=0.
     258          zdtsw1_lev(1)=0.
     259          zdtlw_lev(1)=0.
     260          zdtsw_lev(1)=0.
     261          ztlev(1)=zt(ig,1)
     262
     263          DO l=1,nlayer-1
     264            ! Calculation for the dust storm fraction
     265            zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
     266                         zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
     267                            (pzlay(ig,l+1)-pzlay(ig,l))
     268           
     269            zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
     270                         zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
     271                            (pzlay(ig,l+1)-pzlay(ig,l))
     272            !! Calculation for the background dust fraction
     273            zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
     274                         pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
     275                            (pzlay(ig,l+1)-pzlay(ig,l))
     276           
     277            zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
     278                         pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
     279                            (pzlay(ig,l+1)-pzlay(ig,l))
     280           
     281            ztlev(l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+          &
     282                         zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /      &
     283                            (pzlay(ig,l+1)-pzlay(ig,l))
     284          ENDDO ! DO l=1,nlayer-1
     285        ENDDO ! DO ig=1,ngrid     
     286
     287        !! **********************************************************************
     288        !! 2.3 Calculation of the vertical velocity : extra heating
     289        !! balanced by adiabatic cooling
     290        DO ig=1,ngrid
     291          IF (storm(ig)) THEN       
     292            scheme(ig)=2         
     293            DO l=1,nlayer
     294              deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l))  &
     295                                          -(zdtlw_lev(l)+zdtsw_lev(l))
     296              wrad(ig,l)=-deltahr(ig,l)/(g/cpp+   &
     297                                         max(zdtvert(l),-0.99*g/cpp))       
     298              !! Limit vertical wind in case of lapse rate close to adiabatic
     299              wrad(ig,l)=max(wrad(ig,l),-wmax)
     300              wrad(ig,l)=min(wrad(ig,l),wmax)
     301            ENDDO
     302          ENDIF ! IF (storm(ig))
     303        ENDDO ! DO ig=1,ngrid
     304
     305      ! **********************************************************************
     306      ! 3. Vertical transport
     307      ! **********************************************************************
     308        !! **********************************************************************
     309        !! 3.1 Transport of background dust + storm dust (concentrated)
     310        DO ig=1,ngrid
     311          IF (storm(ig)) THEN
     312            DO l=1,nlayer
     313              mr_dust_mass(ig,l) = zq_dust_mass(ig,l)
     314              mr_dust_number(ig,l) = zq_dust_number(ig,l)
     315              mr_stormdust_mass(ig,l) = zq_dust_mass(ig,l)+ &
     316                                      zq_stormdust_mass(ig,l)/totstormfract(ig)
     317              mr_stormdust_number(ig,l) = zq_dust_number(ig,l)+ &
     318                                      zq_stormdust_number(ig,l)/totstormfract(ig)
     319            ENDDO ! DO l=1,nlayer
     320          ENDIF ! IF (storm(ig))
     321        ENDDO ! DO ig=1,ngrid
     322
     323        !! **********************************************************************
     324        !! 3.2 Vertical transport by a Van Leer scheme
     325        DO ig=1,ngrid
     326          IF (storm(ig)) THEN
     327           
     328            !! Mass of atmosphere in the layer
     329            DO l=1,nlayer
     330               masse(l)=(pplev(ig,l)-pplev(ig,l+1))/g
     331            ENDDO
     332           
     333            !! Mass flux in kg/m2
     334            DO l=1,nlayer
     335               w(l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(l)))*ptimestep
     336            ENDDO
     337           
     338            !! Vector column
     339            DO l=1,nlayer
     340               zq_vl_col(l)= mr_stormdust_mass(ig,l)
     341               zn_vl_col(l)= mr_stormdust_number(ig,l)
     342            ENDDO
     343           
     344            !! Van Leer scheme
     345            wqmass(:)=0.
     346            wqnumber(:)=0.
     347            CALL vl_storm(nlayer,zq_vl_col,2.,   &
     348                  masse,w,wqmass)
     349            CALL vl_storm(nlayer,zn_vl_col,2.,  &
     350                  masse,w,wqnumber)
     351            !! Mass mixing ratio after transport     
     352            mr_stormdust_mass(ig,:) = zq_vl_col(:)
     353            mr_stormdust_number(ig,:) = zn_vl_col(:)
     354                       
     355          ENDIF ! IF storm(ig)
     356        ENDDO ! DO ig=1,ngrid 
     357
     358        !! **********************************************************************
     359        !! 3.3 Re-calculation of the mixing ratios after vertical transport
     360        DO ig=1,ngrid
     361         IF (storm(ig)) THEN
     362           DO l=1,nlayer
     363           
     364             !! General and "healthy" case
     365             IF (mr_stormdust_mass(ig,l).ge.mr_dust_mass(ig,l)) THEN
     366               zq_dust_mass(ig,l) = mr_dust_mass(ig,l)
     367               zq_dust_number(ig,l) = mr_dust_number(ig,l)
     368               zq_stormdust_mass(ig,l) = totstormfract(ig)*(mr_stormdust_mass(ig,l)-mr_dust_mass(ig,l))
     369               zq_stormdust_number(ig,l) = totstormfract(ig)*(mr_stormdust_number(ig,l)-mr_dust_number(ig,l))
     370             !! Particular case: there is less than initial dust mixing ratio after the vertical transport
     371             ELSE
     372               zq_dust_mass(ig,l) = (1.-totstormfract(ig))*mr_dust_mass(ig,l)+totstormfract(ig)*mr_stormdust_mass(ig,l)
     373               zq_dust_number(ig,l) = (1.-totstormfract(ig))*mr_dust_number(ig,l)+totstormfract(ig)*mr_stormdust_number(ig,l)
     374               zq_stormdust_mass(ig,l) = 0.
     375               zq_stormdust_number(ig,l) = 0.
     376             ENDIF
     377             
     378           ENDDO ! DO l=1,nlayer           
     379         ENDIF ! IF storm(ig)
     380        ENDDO ! DO ig=1,ngrid
     381
     382        !! **********************************************************************
     383        !! 3.4 Calculation of the tendencies of the vertical transport
     384        DO ig=1,ngrid
     385         IF (storm(ig)) THEN
     386           DO l=1,nlayer
     387             dqvl_stormdust_mass(ig,l) = (zq_stormdust_mass(ig,l)- &
     388                                   zq(ig,l,igcm_stormdust_mass)) /ptimestep
     389             dqvl_stormdust_number(ig,l) = (zq_stormdust_number(ig,l)- &
     390                                     zq(ig,l,igcm_stormdust_number)) /ptimestep
     391             dqvl_dust_mass(ig,l) = (zq_dust_mass(ig,l)-zq(ig,l,igcm_dust_mass)) /ptimestep
     392             dqvl_dust_number(ig,l) = (zq_dust_number(ig,l)-zq(ig,l,igcm_dust_number)) /ptimestep
     393           ENDDO
     394         ENDIF ! IF storm(ig)
     395        ENDDO ! DO ig=1,ngrid           
     396
     397      ! **********************************************************************
     398      ! 4. Detrainment: stormdust is converted to background dust
     399      ! **********************************************************************
     400        !! **********************************************************************
     401        !! 4.1 Compute the coefficient of detrainmen
     402        DO ig=1,ngrid
     403          DO l=1,nlayer
     404            IF ((max(abs(wrad(ig,l)),abs(wrad(ig,l+1))) .lt.  &
     405                          wmin) .or. (zq_dust_mass(ig,l) .gt.  &
     406                                     10000.*zq_stormdust_mass(ig,l))) THEN
     407               coefdetrain=1.
     408            ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1)))   &
     409                                                       .le. wmax) THEN
     410               coefdetrain=1.*(((1-coefmin)/(wmin-wmax)**2)*     &
     411                   (max(abs(wrad(ig,l)),abs(wrad(ig,l+1)))-wmax)**2 &
     412                                                               +coefmin)
     413            ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1))).gt. wmax ) THEN
     414               coefdetrain=coefmin
     415            ELSE
     416               coefdetrain=coefmin
     417            ENDIF
     418          ENDDO ! DO l=1,nlayer
     419        ENDDO ! DO ig=1,ngrid
    341420       
    342         ! whatever the situation is, we need the vertical velocity computed by
    343         !  the rds scheme
    344         zdtvert(1)=0. !This is the env. lapse rate
    345         DO l=1,nlayer-1
    346           zdtvert(l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
    347           lapserate(ig,l+1)=zdtvert(l+1) !for diagnosing
    348         ENDDO
    349      
    350         ! computing heating rates gradient at boundraies of each layer
    351         ! start from surface
    352         zdtlw1_lev(1)=0.
    353         zdtsw1_lev(1)=0.
    354         zdtlw_lev(1)=0.
    355         zdtsw_lev(1)=0.
    356         ztlev(1)=zt(ig,1)
    357 
    358         DO l=1,nlayer-1
    359           ! Calculation for the dust storm fraction
    360           zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
    361                        zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
    362                           (pzlay(ig,l+1)-pzlay(ig,l))
    363          
    364           zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
    365                        zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
    366                           (pzlay(ig,l+1)-pzlay(ig,l))
    367           !MV18: calculation for the background dust fraction
    368           zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
    369                        pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
    370                           (pzlay(ig,l+1)-pzlay(ig,l))
    371          
    372           zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
    373                        pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
    374                           (pzlay(ig,l+1)-pzlay(ig,l))
    375          
    376           ztlev(l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+          &
    377                        zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /      &
    378                           (pzlay(ig,l+1)-pzlay(ig,l))
    379         ENDDO
    380      
    381         DO l=1,nlayer
    382           deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l))  &
    383                                       -(zdtlw_lev(l)+zdtsw_lev(l))
    384           wadiabatic(ig,l)=-deltahr(ig,l)/(g/cpp+   &
    385                                      max(zdtvert(l),-0.99*g/cpp))
    386      
    387           !limit vertical wind in case of lapse rate close to adiabat
    388           wadiabatic(ig,l)=max(wadiabatic(ig,l),-wlim)
    389           wadiabatic(ig,l)=min(wadiabatic(ig,l),wlim)
    390         ENDDO
     421        !! **********************************************************************
     422        !! 4.2 Calculation of the tendencies of the detrainment
     423        DO ig=1,ngrid
     424          DO l=1,nlayer
     425            dqdet_stormdust_mass(ig,l)=-coefdetrain*zq_stormdust_mass(ig,l) &
     426                                                        /ptimestep
     427            dqdet_stormdust_number(ig,l)=-coefdetrain*zq_stormdust_number(ig,l) &
     428                                                        /ptimestep
     429          ENDDO ! DO l=1,nlayer
     430        ENDDO ! DO ig=1,ngrid
     431           
     432      ! **********************************************************************
     433      ! 5. Final tendencies: vertical transport + detrainment
     434      ! **********************************************************************
     435        DO ig=1,ngrid
     436          DO l=1,nlayer     
     437          pdqrds(ig,l,igcm_stormdust_mass)=dqdet_stormdust_mass(ig,l) &
     438                                                 +dqvl_stormdust_mass(ig,l)
     439          pdqrds(ig,l,igcm_stormdust_number)=dqdet_stormdust_number(ig,l) &
     440                                                 +dqvl_stormdust_number(ig,l)
     441          pdqrds(ig,l,igcm_dust_mass)= -dqdet_stormdust_mass(ig,l) &
     442                                       +dqvl_dust_mass(ig,l)
     443          pdqrds(ig,l,igcm_dust_number)= -dqdet_stormdust_number(ig,l) &
     444                                       +dqvl_dust_number(ig,l)
     445          ENDDO ! DO l=1,nlayer
     446        ENDDO ! DO ig=1,ngrid
     447
     448      ! **********************************************************************
     449      ! 6. To prevent from negative values
     450      ! **********************************************************************
     451        DO ig=1,ngrid
     452          DO l=1,nlayer
     453            IF ((pq(ig,l,igcm_stormdust_mass) &
     454               +pdqrds(ig,l,igcm_stormdust_mass)*ptimestep .le. 0.) .or. &
     455              (pq(ig,l,igcm_stormdust_number) &
     456               +pdqrds(ig,l,igcm_stormdust_number)*ptimestep .le. 0.)) THEN
     457               pdqrds(ig,l,igcm_stormdust_mass)=-pq(ig,l,igcm_stormdust_mass)/ptimestep
     458               pdqrds(ig,l,igcm_stormdust_number)=-pq(ig,l,igcm_stormdust_number)/ptimestep
     459            ENDIF
     460          ENDDO ! nlayer
     461        ENDDO ! DO ig=1,ngrid
     462
     463        DO ig=1,ngrid
     464          DO l=1,nlayer           
     465            IF ((pq(ig,l,igcm_dust_mass) &
     466               +pdqrds(ig,l,igcm_dust_mass)*ptimestep .le. 0.) .or. &
     467              (pq(ig,l,igcm_dust_number) &
     468               +pdqrds(ig,l,igcm_dust_number)*ptimestep .le. 0.)) THEN
     469               pdqrds(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)/ptimestep
     470               pdqrds(ig,l,igcm_dust_number)=-pq(ig,l,igcm_dust_number)/ptimestep
     471            ENDIF
     472          ENDDO ! nlayer
     473        ENDDO ! DO ig=1,ngrid
    391474       
    392           !! **********************************************************************
    393           !! 2.2 case 2: daytime slope wind scheme     
    394           IF ((mu0(ig) .gt. mu0lim) .and. .not. storm(ig)) then
    395             scheme(ig)=2
    396             alpha(ig) = alpha_hmons(ig)
    397             ! interpolate the env. temperature above the mountain top
    398             call intep_vtemp(nlayer,hmons(ig),zt(ig,:),pzlay(ig,:), &
    399                                envtemp,slpbg(ig))
    400             envt(ig,:)=envtemp(:) !for diagnosis
    401  
    402             ! second: estimate the vertical velocity at boundraies of each layer
    403             !wspeed(ig,1)=0. ! at surface, already initialized
    404  
    405             !!!!!!!the first layer of atmosphere!!!!!!!!!!
    406             IF (slpbg(ig) .gt. 0.) THEN !only postive buoyancy
    407                ! if slpbg(ig) lt 0, means the slope flow is colder than env. (night or
    408                ! early morning ?)  !!!!!!!!!!!
    409                ! ideal method
    410               !w1(ig)=-sqrt(g*slpbg(ig)/zt(ig,1)*hmons(ig))*sin(zsig(ig))
    411                ! new scheme, simply proportional to temperature and
    412                ! mountain height
    413                w0(ig)=-1.5e-4*g*slpbg(ig)/zt(ig,1)*hmons(ig)
    414                ! otherwise, w0(ig) =0.
    415                wspeed(ig,2)=w0(ig)
    416             ELSE
    417                wspeed(ig,2)=wadiabatic(ig,2) !! for early morning ?
    418             ENDIF
    419        
    420             ! prepare the integration, NOTE: if w is too small, may have artificials
    421             IF (abs(wspeed(ig,2)) .lt. 0.01 ) &
    422                                  wspeed(ig,2)=sign(0.01,wspeed(ig,2))
    423  
    424             newzt(ig,1)= zt(ig,1) !temperature of the first layer atmosphere
    425                                   !  above the mountain top (radiative
    426                                   !  equilibrium on Mars)
    427 
    428             ! estimate the vertical velocities
    429             ! if w0(ig) >= 0, means downward motion, no upslope winds, we switch to
    430             ! assume that the extra heating integrally convert to
    431             ! vertical motion.
    432             if ( w0(ig) .ge. 0 ) then  !! normal, it is impossible,
    433                                        !!  because mu(ig)>0.1 here
    434                do l=3,nlayer
    435                   wspeed(ig,l)=wadiabatic(ig,l)
    436                enddo
    437             else
    438             ! estimate the velocities by taking into account the heating due
    439             ! to storm dust, the cooling due to vertical motion ...
    440             !!!!!!!!!!!the simple scheme!!!!!!!!!
    441                do l=2,nlayer-1
    442                  !if superadiabatic layer
    443                  if ( zdtvert(l) .lt. -g/cpp) then
    444                     !case 1
    445                     ! test, also decrease adiabatically ?
    446                    !newzt(ig,l)= &
    447                    !        zt(ig,l-1)-g/cpp*(pzlay(ig,l)-pzlay(ig,l-1))
    448                     newzt(ig,l)=zt(ig,l)
    449                    !wspeed(ig,l+1)=wspeed(ig,l)
    450                  else
    451                     !not superadiabatic
    452                     newzt(ig,l)=newzt(ig,l-1)+(deltahr(ig,l)/ &
    453                                   (-wspeed(ig,l))-g/cpp)*        &
    454                                   (pzlay(ig,l)-pzlay(ig,l-1))
    455                  endif ! end of if superadiabatic or not
    456  
    457                 !wtemp(ig,l+1)=wspeed(ig,l)**2+2.*g*(pzlev(ig,l+1) &
    458                 !              -pzlev(ig,l))*(k1*(newzt(ig,l)       &
    459                 !              -envtemp(l))/envtemp(l))
    460                  wtemp(ig,l+1)=(1.-2.*k1*(pzlev(ig,l+1)-pzlev(ig,l)))*&
    461                                 wspeed(ig,l)**2+2.*k2*g*(pzlev(ig,l+1) &
    462                                -pzlev(ig,l))*((newzt(ig,l)       &
    463                                -envtemp(l))/envtemp(l))
    464  
    465                  if (wtemp(ig,l+1) .gt. 0.) then
    466                    !case 2
    467                    wspeed(ig,l+1)=-sqrt(wtemp(ig,l+1))
    468 
    469                    ! if |wspeed| < |wadiabatic| then go to wadiabatic
    470                    if (wspeed(ig,l+1) .gt. wadiabatic(ig,l+1)) then
    471                      do ll=l,nlayer-1
    472                        newzt(ig,ll)=envtemp(ll)
    473                        wspeed(ig,ll+1)=wadiabatic(ig,ll+1)
    474                      enddo
    475                      exit
    476                    endif
    477 
    478                    ! avoid artificials
    479                    if (abs(wspeed(ig,l+1)) .lt. 0.01 ) &
    480                               wspeed(ig,l+1)=sign(0.01,wspeed(ig,l+1))
    481        
    482                  else if (l .lt. nlayer) then
    483                    !case 3
    484                    do ll=l,nlayer-1
    485                      newzt(ig,ll)=envtemp(ll)
    486                      wspeed(ig,ll+1)=wadiabatic(ig,ll+1)
    487                    enddo !overshot
    488                    exit
    489          
    490                  else
    491                    wspeed(ig,l+1)=0.
    492                  endif
    493        
    494                enddo
    495        
    496             endif !w0
    497  
    498          ELSE
    499           !! **********************************************************************
    500           !! 2.3 case 3: storm=true
    501           if (storm(ig)) then
    502               scheme(ig)=3
    503               alpha(ig) = totstormfract(ig)
    504               do l=1,nlayer
    505                   wspeed(ig,l)=wadiabatic(ig,l)
    506               enddo
    507            endif ! storm=1
    508  
    509          ENDIF ! rds or slp
    510 
    511 
    512         !!!!!!!! estimate the amount of dust for diagnostics
    513         DO l=1,nlayer
    514             ! transfer background dust + storm dust (concentrated)
    515             zqstormdustslp(ig,l) =zqi_mass(ig,l)+ &
    516                                       zqstorm_mass(ig,l)/alpha(ig)
    517             znstormdustslp(ig,l) =zqi_number(ig,l)+ &
    518                                     zqstorm_number(ig,l)/alpha(ig)
    519             zqdustslp(ig,l) = zqi_mass(ig,l)
    520             zndustslp(ig,l) = zqi_number(ig,l)
    521      
    522             rdsdustqvl0(ig,l)=zqstormdustslp(ig,l) !for diagnosis
    523         ENDDO
    524 
    525       ! **********************************************************************
    526       ! 3. Vertical transport
    527       ! **********************************************************************
    528         do l=1,nlayer
    529              masse(l)=(pplev(ig,l)-pplev(ig,l+1))/g
    530         enddo
    531         !Estimation of "dtmax" (s) to be used for vertical transport   
    532         dtmax=ptimestep
    533         !secu is a margin on subtimstep to avoid dust crossing many layers
    534         do l=2,nlayer
    535           if (wspeed(ig,l).lt.0.) then       ! case up
    536              dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/  &
    537                                (secu*abs(wspeed(ig,l))))
    538           else if (wspeed(ig,l).gt.0.) then
    539              dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/  &
    540                                (secu*abs(wspeed(ig,l))))
    541           endif
    542         enddo
    543 
    544         nsubtimestep= int(ptimestep/dtmax)
    545         subtimestep=ptimestep/float(nsubtimestep)
    546        
    547         do l=1,nlayer
    548            wrds(l)=wspeed(ig,l)*pplev(ig,l)/(r*ztlev(l)) &
    549                     *subtimestep
    550         enddo
    551      
    552         do l=1,nlayer
    553            zqstorm_mass_col(l)= zqstormdustslp(ig,l) !zqstorm_mass(ig,l)
    554            zqstorm_number_col(l)=znstormdustslp(ig,l) ! zqstorm_number(ig,l)
    555         enddo
    556        
    557         do tsub=1,nsubtimestep
    558            wqrdsmass(:)=0.
    559            wqrdsnumber(:)=0.
    560            CALL vl_storm(nlayer,zqstorm_mass_col,2.,   &
    561                  masse,wrds ,wqrdsmass)
    562            CALL vl_storm(nlayer,zqstorm_number_col,2.,  &
    563                  masse,wrds ,wqrdsnumber)
    564         enddo
    565      
    566         !!!!!generate the "extra" dust
    567         do l=1,nlayer
    568            rdsdustqvl1(ig,l)=zqstorm_mass_col(l) ! for diagnosis
    569 
    570           ! extra dust = storm dust
    571           !zqdustslp(ig,l)=zqi_mass(ig,l) !(1.-alpha(ig))*zqi_mass(ig,l)
    572           !zndustslp(ig,l)=zqi_number(ig,l) !(1.-alpha(ig))*zqi_number(ig,l)
    573           !zqstorm_mass_col(l)=alpha(ig)*zqstorm_mass_col(l)
    574           !zqstorm_number_col(l)=alpha(ig)*zqstorm_number_col(l)
    575 
    576            !with compensation
    577            if (zqstorm_mass_col(l) .lt. zqi_mass(ig,l)  ) then
    578               ! the following two equations are easier to understand
    579               zqdustslp(ig,l)=(1.-alpha(ig))*zqi_mass(ig,l)+alpha(ig)* &
    580                          zqstorm_mass_col(l)
    581               zndustslp(ig,l)=(1.-alpha(ig))*zqi_number(ig,l)+alpha(ig)&
    582                          *zqstorm_number_col(l)
    583               !with a bug, should be  zqi+alpha****
    584               !zqdustslp(ig,l)=zqi_mass(ig,l)-alpha(ig)* &
    585               !           (zqstorm_mass_col(l)-zqi_mass(ig,l))
    586               !zndustslp(ig,l)=zqi_number(ig,l)-alpha(ig)* &
    587               !           (zqstorm_number_col(l)-zqi_number(ig,l))
    588               zqstorm_mass_col(l)=0. 
    589               zqstorm_number_col(l)=0.
    590            else
    591               zqstorm_mass_col(l)=alpha(ig)* &
    592                          (zqstorm_mass_col(l)-zqi_mass(ig,l))
    593               zqstorm_number_col(l)=alpha(ig)* &
    594                          (zqstorm_number_col(l)-zqi_number(ig,l))
    595               ! the mass mixing ratio of environmental dust doesn't change.
    596            endif
    597            !diagnose
    598            stormdust_m1(ig,l)=zqstorm_mass_col(l)
    599         enddo
    600 
    601         !=======================================================================
    602         ! calculate the tendencies due to vertical transport
    603         do l=1,nlayer
    604            ! tendencies due to vertical transport
    605            zdqvlstorm_mass(ig,l)= (zqstorm_mass_col(l)- &
    606                                    zqstorm_mass(ig,l)) /ptimestep
    607            zdqvlstorm_number(ig,l)= (zqstorm_number_col(l)- &
    608                                  zqstorm_number(ig,l)) /ptimestep
    609 
    610            zdqenv_mass(ig,l)=(zqdustslp(ig,l)-zqi_mass(ig,l))/ptimestep
    611            zdqenv_number(ig,l)=(zndustslp(ig,l)-zqi_number(ig,l))  &
    612                                 /ptimestep
    613 
    614            ! chao for output only
    615           !qstormdustvl1(ig,l)=zqstorm_mass_col(l)
    616           !nstormdustvl1(ig,l)=zqstorm_number_col(l)
    617           !stormdust_m_col1(ig)=stormdust_m_col1(ig)+zqstorm_mass_col(l)  &
    618           !                      *(pplev(ig,l)-pplev(ig,l+1))/g
    619           !rdsdustqvl1(ig,l)=zqstorm_mass_col(l)
    620         enddo
    621 
    622       ! **********************************************************************
    623       ! 4. Detrainment: convert dust storm to background dust
    624       ! **********************************************************************
    625         do l=1,nlayer
    626              ! compute the coefficient of detrainment
    627           if ((max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))) .lt.  &
    628                         detrainlim) .or. (zqdustslp(ig,l) .gt.  &
    629                                    10000.*zqstorm_mass_col(l))) then
    630              coefdetrain=1.
    631           else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))   &
    632                                                      .le. wlim) then
    633                   ! case where detrainment depends on vertical wind
    634            ! coefdetrain=0.5*(((1-coefmin)/(detrainlim-3.)**2)*     &
    635            !     (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-3.)**2 &
    636            !                                               +coefmin)
    637              coefdetrain=1.*(((1-coefmin)/(detrainlim-wlim)**2)*     &
    638                  (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-wlim)**2 &
    639                                                              +coefmin)
    640            !coefdetrain=0.5
    641           else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))).gt. 10. )&
    642                   then
    643              coefdetrain=0.025
    644           else
    645              coefdetrain=coefmin
    646           endif
    647 
    648           detrainment(ig,l)=coefdetrain !for diagnose
    649 
    650           ! Calculate tendancies corresponding to pdq after detrainement
    651           ! pdqdet = tendancy corresponding to detrainment only
    652           zdqdetstorm_mass(ig,l)=-coefdetrain*zqstorm_mass_col(l) &
    653                                                         /ptimestep
    654           zdqdetstorm_number(ig,l)=-coefdetrain*zqstorm_number_col(l) &
    655                                                         /ptimestep
    656      
    657           ! pdqrds ( tendancy corresponding to vertical transport and
    658           ! detrainment) = zdqvlstorm + pdqdet
    659           pdqrds(ig,l,igcm_stormdust_mass)=zdqdetstorm_mass(ig,l) &
    660                                                  +zdqvlstorm_mass(ig,l)
    661           pdqrds(ig,l,igcm_stormdust_number)=zdqdetstorm_number(ig,l) &
    662                                                +zdqvlstorm_number(ig,l)
    663           pdqrds(ig,l,igcm_dust_mass)= zdqenv_mass(ig,l) &
    664                              -zdqdetstorm_mass(ig,l)
    665           pdqrds(ig,l,igcm_dust_number)= zdqenv_number(ig,l) &
    666                              -zdqdetstorm_number(ig,l)
    667 
    668           !diagnose
    669           stormdust_m2(ig,l)=zqstorm_mass_col(l)-coefdetrain*zqstorm_mass_col(l)
    670         enddo ! nlayer
    671 !       endif
    672       !=======================================================================
    673       enddo ! end do ig=1,ngrid
    674  
    675 !     !chao check conservation here
    676 !     do l=1,nlayer
    677 !       do ig=1,ngrid
    678 !         totdust0(ig)=totdust0(ig)+                &
    679 !            zq(ig,l,igcm_stormdust_mass)*        &
    680 !            ((pplev(ig,l) - pplev(ig,l+1)) / g)  &
    681 !            +  zq(ig,l,igcm_dust_mass)*          &
    682 !            ((pplev(ig,l) - pplev(ig,l+1)) / g)
    683 
    684 !         totdust1(ig)=totdust1(ig)+                &
    685 !            (zq(ig,l,igcm_stormdust_mass) + &
    686 !             pdqrds(ig,l,igcm_stormdust_mass)*ptimestep)*        &
    687 !            ((pplev(ig,l) - pplev(ig,l+1)) / g)  &
    688 !            +  ( zq(ig,l,igcm_dust_mass)+          &
    689 !             pdqrds(ig,l,igcm_dust_mass)*ptimestep)*        &
    690 !            ((pplev(ig,l) - pplev(ig,l+1)) / g)
    691 !       enddo
    692 !     enddo
    693      
    694 !       call writediagfi(ngrid,'totdust0','total dust before rds', &
    695 !              ' ',2,totdust0)
    696 !       call writediagfi(ngrid,'totdust1','total dust after rds', &
    697 !              ' ',2,totdust1)
    698         !output for diagnosis
    699         call WRITEDIAGFI(ngrid,'detrainment',         &
    700                         'coefficient of detrainment',' ',3,detrainment)
    701        !call WRITEDIAGFI(ngrid,'qstormvl1','mmr of stormdust after rds_vl', &
    702        !   &                        'kg/kg',3,qstormdustvl1)
     475!=======================================================================
     476! WRITEDIAGFI
     477
    703478        call WRITEDIAGFI(ngrid,'lapserate','lapse rate in the storm', &
    704479           &                        'k/m',3,lapserate)
    705480        call WRITEDIAGFI(ngrid,'deltahr','extra heating rates', &
    706481           &                        'K/s',3,deltahr)
    707         call WRITEDIAGFI(ngrid,'wold', &
    708                              'wind generated from adiabatic colling ', &
    709            &                        'm/s',3,wadiabatic)
    710         call WRITEDIAGFI(ngrid,'newzt','perturbated temperature', &
    711            &                        'K/s',3,newzt)
    712         call WRITEDIAGFI(ngrid,'zt','unperturbated temperature', &
    713            &                        'K/s',3,zt)
    714         call WRITEDIAGFI(ngrid,'wtemp','under square root', &
    715            &                        'K/s',3,wtemp)
    716        !call WRITEDIAGFI(ngrid,'stormdust_m_col1','mass of stormdust after rds_vl', &
    717        !   &                        'kg/kg',2,stormdust_m_col1)
    718        !call WRITEDIAGFI(ngrid,'temprds','temp for calculating zdtvert', &
    719        !   &                        'k',3,temprds)
    720         call WRITEDIAGFI(ngrid,'stormdust_m0','mass col  of stormdust before rds_vl', &
    721            &                        'kg/kg',3,stormdust_m0)
    722         call WRITEDIAGFI(ngrid,'stormdust_m1','mass col  of stormdust after rds_vl', &
    723            &                        'kg/kg',3,stormdust_m1)
    724         call WRITEDIAGFI(ngrid,'stormdust_m2','mass col  of stormdust after rds_vl', &
    725            &                        'kg/kg',3,stormdust_m2)
    726      
    727       ! call writediagfi(ngrid,'wslp','estimated slope winds','m/s',3,wslp)
    728       ! call writediagfi(ngrid,'wslp2','estimated slope winds2','m/s',3,wslp2)
    729       ! call writediagfi(ngrid,'zhb','estimated slope winds2','m/s',3,zhb)
    730       ! call writediagfi(ngrid,'bruntf','bouyancy frequency',' ',3,bruntf)
    731       ! call writediagfi(ngrid,'slpdepth','slope depth','m',2,slpdepth)
    732       ! call writediagfi(ngrid,'slpu','perbulation u','m/s',3,slpu)
    733       ! call writediagfi(ngrid,'slpzh','perbulation zh',' ',3,slpzh)
    734       ! call writediagfi(ngrid,'zqslp','zq in rocketduststorm','ikg/kg',3, &
    735       !                     zq(:,:,igcm_dust_mass))
    736       ! call writediagfi(ngrid,'zrdsqslp','zq in rocketduststorm','ikg/kg',3, &
    737       !                     zq(:,:,igcm_stormdust_mass))
    738       ! call writediagfi(ngrid,'wslplev','estimated slope winds','m/s',3,wslplev)
    739       ! call writediagfi(ngrid,'slope','identified slope wind effect',' ',2,slope)
    740         call writediagfi(ngrid,'w0','max of slope wind',' ',2,w0)
    741 !       call writediagfi(ngrid,'w1','max of slope wind',' ',2,w1)
    742         call writediagfi(ngrid,'mu0','cosine of solar incidence angle',&
    743                                                            ' ',2,mu0)
    744 !       call writediagfi(ngrid,'storm','identified rocket dust storm',&
    745 !                                                  ' ',2,real(storm))
    746482        call writediagfi(ngrid,'scheme','which scheme',&
    747483                                                   ' ',2,real(scheme))
    748         call writediagfi(ngrid,'alpha','coefficient alpha',' ',2,alpha)
    749       ! call writediagfi(ngrid,'q2rds','alpha zq',' ', &
    750       !                         3,q2rds)
    751         call writediagfi(ngrid,'rdsdustqvl0','not vl storm slp',       &
    752                                               'kg/kg',3,zqstormdustslp)
    753         call writediagfi(ngrid,'rdsdustqvl1','vled storm slp',         &
    754                                                  'kg/kg',3,rdsdustqvl1)
    755         call writediagfi(ngrid,'dustqvl0','not vl  slp',               &
    756                                                     'kg/kg',3,zqi_mass)
    757         call writediagfi(ngrid,'dustqvl1','vled  slp',                 &
    758                                                    'kg/kg',3,zqdustslp)
    759       ! call WRITEDIAGFI(ngrid,'lmax_th2', &
    760       !                  'hauteur du thermique','point', &
    761       !                             2,real(lmax_th(:)))
    762       ! call WRITEDIAGFI(ngrid,'zmax_th2', &
    763       !                  'hauteur du thermique','m', &
    764       !                             2,zmax_th)
    765       ! call writediagfi(ngrid,'lslope','lenght of slope',' ',2,lslope)
    766       ! call writediagfi(ngrid,'hmon','identified slope wind effect',' ',2,hmon)
    767         call writediagfi(ngrid,'envt','interpolated env. temp.',       &
    768                                                            'K',3,envt)
    769         call writediagfi(ngrid,'hmons','identified slope wind effect', &
    770                                                            ' ',2,hmons)
    771       ! call writediagfi(ngrid,'slpbg','temp. diff along slope',       &
    772       !                                                    ' ',2,slpbg)
    773       ! call writediagfi(ngrid,'zmea','identified slope wind effect',  &
    774       !                                                    ' ',2,zmea)
    775       ! call writediagfi(ngrid,'zsig','identified slope wind effect',  &
    776       !                                                    ' ',2,zsig)
    777       ! call writediagfi(ngrid,'zhslpenv','difference of zh above mons',' ',2,zhslpenv)
    778       ! call writediagfi(ngrid,'lref','identified slope wind effect',' ',2,real(lref))
     484
    779485        END SUBROUTINE rocketduststorm
    780486
    781 !********************************************************************************
    782 !********************************************************************************
     487!=======================================================================
     488! VAN LEER
     489
    783490      SUBROUTINE vl_storm(nlay,q,pente_max,masse,w,wq)
    784491!
    … …  
    801508!   Arguments:
    802509!   ----------
    803       integer,intent(in) :: nlay ! number of atmospheric layers
    804       real masse(nlay),pente_max
     510      INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers
     511      REAL masse(nlay),pente_max
    805512      REAL q(nlay)
    806513      REAL w(nlay)
    … …  
    813520!
    814521
    815       real dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
    816       real newmasse
    817       real sigw, Mtot, MQtot
    818       integer m
     522      REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
     523      REAL newmasse
     524      REAL sigw, Mtot, MQtot
     525      INTEGER m
    819526
    820527      REAL      SSUM,CVMGP,CVMGT
    821       integer ismax,ismin
     528      INTEGER ismax,ismin
    822529
    823530
    … …  
    847554      dzq(1)=0.
    848555      dzq(nlay)=0.
    849        
    850 !      write(*,*),'TB14 wq before up',wq(1,:)
    851 !      write(*,*),'TB14 q before up',q(1,:)
     556
    852557! ---------------------------------------------------------------
    853558!   .... calcul des termes d'advection verticale  .......
    … …  
    891596             if (m.gt.0) then
    892597                sigw=(w(l+1)+Mtot)/masse(m)
    893                 wq(l+1)= (MQtot + (-w(l+1)-Mtot)*         &
     598                wq(l+1)= -(MQtot + (-w(l+1)-Mtot)*         &
    894599                (q(m)-0.5*(1.+sigw)*dzq(m))  )
    895600             else
    896 ! new
    897601                w(l+1) = -Mtot
    898602                wq(l+1) = -MQtot
    899 !                write(*,*) 'TB14 MQtot = ',MQtot,l
    900 
    901 ! old
    902 !                wq(l+1)= (MQtot + (-w(l+1)-Mtot)*qm(1))
    903 !                write(*,*) 'a rather weird situation in vlz_fi !'
    904 !                stop
    905603             end if
    906604          endif
    … …  
    913611
    914612       enddo
    915 
    916 !       write(*,*),'TB14 masse',masse(1,:)
    917 !       write(*,*),'TB14 wq before down after up',wq(1,:)
    918 !       write(*,*),'TB14 q before down',q(1,:)
    919613     
    920614!      2) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION)     
    … …  
    967661
    968662       enddo
    969       end subroutine vl_storm
    970 !********************************************************************************
    971       SUBROUTINE intep_vtemp(nlayer,hm,temp,zlay,envtemp,slpb)
    972 
    973        USE comcstfi_h, only: g,cpp
    974 
     663       
     664      END SUBROUTINE vl_storm
     665
     666!=======================================================================
     667! Initialization of the module variables
     668
     669       subroutine ini_rocketduststorm_mod(ngrid)
     670       
    975671       implicit none
    976        ! this subroutine is using for obtaining the environmental
    977        ! temperature profile when a subgrid mountain exists.
    978 
    979 
    980 !      input:
    981        integer,intent(in) :: nlayer
    982        real,intent(in) :: hm  ! the height of mountain
    983        real,intent(in) :: temp(nlayer)   !large scale temp. profile of one mesh
    984        real,intent(in) :: zlay(nlayer)   ! height at the middle of each layer
    985 
    986 !      output:
    987        real,intent(out) :: envtemp(nlayer) ! environment temperature
    988        real,intent(out) :: slpb !the temperature difference between slope and
    989                                 ! env. at the half height of mountain
    990 
    991 !      local variables
    992        integer l,il,tmpl
    993        integer lnew  !the layer of atmosphere above the mountain
    994                              ! corresponding to the env. (for buoyancy
    995                              ! calc. )
    996        real newh(nlayer)     !height at the middle of each layer
    997                              ! account for the exist of mountain
    998 !      real g,cpp
    999        real halfh !  half the height of a mountain
    1000 
    1001        !initilize the array
    1002        lnew=0
    1003        newh(:)=0.
    1004        envtemp(:)=0.
    1005        tmpl=1 
    1006 
    1007        do l=1,nlayer
    1008          newh(l)=hm+zlay(l)
    1009          
    1010          do il=tmpl,nlayer-1 !MV18: added the -1
    1011             if (newh(l) .ge. zlay(il) .and. newh(l) .lt. zlay(il+1))then
    1012               ! find the corresponding layer
    1013               lnew=il
    1014 
    1015               ! interpolate
    1016               envtemp(l) = temp(il)+(newh(l)-zlay(lnew))*&
    1017                           (temp(il+1)-temp(il))/(zlay(il+1)-zlay(il))
    1018 
    1019               exit !go to the next layer
    1020             else if (newh(l) .ge. zlay(nlayer)) then
    1021               ! higher than the highest layer
    1022               lnew=nlayer
    1023               envtemp(l)=temp(nlayer)   
    1024 
    1025             endif
    1026          enddo
    1027          ! this can accelerate the loop
    1028          tmpl=lnew
    1029 
    1030        enddo
    1031 
    1032        halfh=0.5*hm
    1033        if (halfh .le. zlay(1) ) then
    1034          slpb=0.
    1035        else if (halfh .gt. zlay(nlayer)) then
    1036          !normally, impossible for halfh gt zlay(l), anyway...
    1037          tmpl=nlayer
    1038          !difference between surface and atmosphere (env.)
    1039          slpb=temp(1)-(temp(nlayer-1)+((halfh-zlay(nlayer-1))* &
    1040               (temp(tmpl)-temp(tmpl-1))/(zlay(tmpl)-zlay(tmpl-1))))
    1041        else
    1042          do l=1,nlayer-1
    1043            if ((halfh .gt. zlay(l)) .and. (halfh .le. zlay(l+1)))then
    1044              tmpl= l
    1045              slpb=temp(1)-(temp(tmpl)+(halfh-zlay(tmpl))*  &
    1046                  (temp(tmpl+1)-temp(tmpl))/(zlay(tmpl+1)-zlay(tmpl)))
    1047            endif
    1048          enddo
    1049        endif         
    1050        end subroutine intep_vtemp
    1051 
    1052        ! initialization module variables
    1053        subroutine ini_rocketduststorm_mod(ngrid)
     672       
     673       integer, intent(in) :: ngrid
     674       
     675       allocate(dustliftday(ngrid))
     676       
     677       end subroutine ini_rocketduststorm_mod
     678       
     679       subroutine end_rocketduststorm_mod
    1054680       
    1055681       implicit none
    1056682       
    1057        integer, intent(in) :: ngrid
    1058        
    1059        allocate(dustliftday(ngrid))
    1060        allocate(alpha_hmons(ngrid))
    1061        
    1062        end subroutine ini_rocketduststorm_mod
    1063        
    1064        subroutine end_rocketduststorm_mod
    1065        
    1066        implicit none
    1067        
    1068683       if (allocated(dustliftday)) deallocate(dustliftday)
    1069        if (allocated(alpha_hmons)) deallocate(alpha_hmons)
    1070684
    1071685       end subroutine end_rocketduststorm_mod       
Note: See TracChangeset for help on using the changeset viewer.