Changeset 3292


Ignore:
Timestamp:
Apr 5, 2024, 8:34:42 AM (9 months ago)
Author:
emillour
Message:

Mars PCM:
Code cleanup to prepare the addition of new schemes for dust lifting by
wind stress and dust devils. Renamed "dustlift" routine "dust_windstress_lift"
and made it a module; also made dustdevil a module.
EM

Location:
trunk/LMDZ.MARS
Files:
3 edited
2 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3290 r3292  
    45944594== 21/03/2024 == Jliu
    45954595Found and fixed another small bug in photochemistry related to commit r3289.
     4596
     4597== 05/04/2024 == EM
     4598Code cleanup to prepare the addition of new schemes for dust lifting by
     4599wind stress and dust devils. Renamed "dustlift" routine "dust_windstress_lift"
     4600and made it a module; also made dustdevil a module.

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

    r3291 r3292  
    1       SUBROUTINE dustlift(ngrid,nlay,nq,rho,
    2      $                  pcdh_true,pcdh,co2ice,
    3      $                  dqslift)
     1MODULE dust_windstress_lift_mod
     2
     3IMPLICIT NONE
     4
     5INTEGER,SAVE :: windstress_lift_scheme ! wind stress lifting scheme number
     6! 0 (default): old scheme
     7
     8!$OMP THREADPRIVATE(windstress_lift_scheme)
     9
     10CONTAINS
     11
     12  SUBROUTINE dust_windstress_lift(ngrid,nlay,nq,rho, &
     13                                pcdh_true,pcdh,co2ice, &
     14                                dqslift)
    415
    516#ifndef MESOSCALE
    6       use tracer_mod, only: alpha_lift, radius
     17  use tracer_mod, only: alpha_lift, radius
    718#else
    8       use tracer_mod, only: alpha_lift, radius,
    9      &                      igcm_dust_mass, igcm_dust_number,
    10      &                      ref_r0,r3n_q
     19  use tracer_mod, only: alpha_lift, radius, &
     20                      igcm_dust_mass, igcm_dust_number, &
     21                      ref_r0,r3n_q
    1122#endif
    12       USE comcstfi_h
    13       IMPLICIT NONE
     23  use comcstfi_h, only: g
     24  USE ioipsl_getin_p_mod, ONLY : getin_p
     25  IMPLICIT NONE
    1426
    15 c=======================================================================
    16 c
    17 c  Dust lifting by surface winds
    18 c  Computing flux to the middle of the first layer
    19 c  (Called by vdifc)
    20 c
    21 c=======================================================================
     27!=======================================================================
     28!
     29!  Dust lifting by surface winds
     30!  Computing flux to the middle of the first layer
     31!  (Called by vdifc)
     32!
     33!=======================================================================
    2234
    23 c-----------------------------------------------------------------------
    24 c   declarations:
    25 c   -------------
     35!
     36!   arguments:
     37!   ----------
    2638
    27 c
    28 c   arguments:
    29 c   ----------
     39  integer,intent(in) :: ngrid, nlay, nq 
     40  real,intent(in) ::  rho(ngrid)  ! density (kg.m-3) at surface
     41  real,intent(in) ::  pcdh_true(ngrid) ! Cd
     42  real,intent(in) ::  pcdh(ngrid) ! Cd * |V|
     43  real,intent(in) ::  co2ice(ngrid) ! surface CO2 ice (kg/m2)
    3044
    31 c   INPUT
    32       integer ngrid, nlay, nq 
    33       real rho(ngrid)  ! density (kg.m-3) at surface
    34       real pcdh_true(ngrid) ! Cd
    35       real pcdh(ngrid) ! Cd * |V|
    36       real co2ice(ngrid)
     45  real,intent(out) ::  dqslift(ngrid,nq) !surface dust flux to mid-layer (<0 when lifing)
    3746
    38 c   OUTPUT
    39       real dqslift(ngrid,nq) !surface dust flux to mid-layer (<0 when lifing)
    40 c     real pb(ngrid,nlay) ! diffusion to surface coeff.
     47!   local:
     48!   ------
     49  INTEGER :: ig,iq
     50  REAL :: fhoriz(ngrid)  ! Horizontal dust flux
     51  REAL :: ust,us
     52  REAL,SAVE :: stress_seuil=0.0225 ! stress lifting threshold (N.m2)
     53!$OMP THREADPRIVATE(stress_seuil)
    4154
    42 c   local:
    43 c   ------
    44       INTEGER ig,iq
    45       REAL fhoriz(ngrid)  ! Horizontal dust flux
    46       REAL ust,us
    47       REAL stress_seuil
    48       SAVE stress_seuil
    49       DATA stress_seuil/0.0225/   ! stress seuil soulevement (N.m2)
    50      
    51 !$OMP THREADPRIVATE(stress_seuil)
     55  LOGICAL,SAVE :: firstcall=.true.
     56!$OMP THREADPRIVATE(firstcall)
     57
     58  character(len=20),parameter :: rname="dust_windstress_lift"
    5259
    5360#ifdef MESOSCALE
    … …  
    5562!!!! AS: ... stress for lifting
    5663!!!! AS: you have to compile with -DMESOSCALE to do so
    57       REAL alpha
    58       REAL r0_lift
    59       INTEGER ierr
    60       REAL ulim
    61         OPEN(99,file='stress.def',status='old',form='formatted'
    62      .   ,iostat=ierr)
    63         !!! no file => default values
    64         IF(ierr.EQ.0) THEN
    65           READ(99,*) ulim !ulim = sqrt(stress_seuil/rho) avec rho = 0.02.
    66                           !prendre ulim = 1.061 m/s pour avoir stress_seuil = 0.0225
    67           READ(99,*) alpha
    68           stress_seuil = 0.02 * ulim * ulim
    69           write(*,*) 'USER-DEFINED threshold: ', stress_seuil, alpha
    70           CLOSE(99)
    71           alpha_lift(igcm_dust_mass) = alpha
    72           r0_lift = radius(igcm_dust_mass) / ref_r0
    73           alpha_lift(igcm_dust_number)=r3n_q*
    74      &                        alpha_lift(igcm_dust_mass)/r0_lift**3
    75           write(*,*) 'set dust number: ', alpha_lift(igcm_dust_number)
    76         ENDIF
     64  REAL :: alpha
     65  REAL :: r0_lift
     66  INTEGER :: ierr
     67  REAL :: ulim
     68
     69    OPEN(99,file='stress.def',status='old',form='formatted',iostat=ierr)
     70    !!! no file => default values
     71    IF(ierr.EQ.0) THEN
     72      READ(99,*) ulim !ulim = sqrt(stress_seuil/rho) avec rho = 0.02.
     73                      !prendre ulim = 1.061 m/s pour avoir stress_seuil = 0.0225
     74      READ(99,*) alpha
     75      stress_seuil = 0.02 * ulim * ulim
     76      write(*,*) 'USER-DEFINED threshold: ', stress_seuil, alpha
     77      CLOSE(99)
     78      alpha_lift(igcm_dust_mass) = alpha
     79      r0_lift = radius(igcm_dust_mass) / ref_r0
     80      alpha_lift(igcm_dust_number)=r3n_q* &
     81                              alpha_lift(igcm_dust_mass)/r0_lift**3
     82      write(*,*) 'set dust number: ', alpha_lift(igcm_dust_number)
     83    ENDIF
    7784#endif
    7885
    79 c     ---------------------------------
    80 c     Computing horizontal flux: fhoriz
    81 c     ---------------------------------
     86    if (firstcall) then
     87      ! get wind stress lifting scheme number (default 0)
     88      windstress_lift_scheme=0 ! default
     89      call getin_p("windstress_lift_scheme",windstress_lift_scheme)
     90     
     91      ! sanity check on available windstress_lift_scheme values
     92      if (windstress_lift_scheme.ne.0) then
     93        write(*,*) trim(rname)//" wrong value for windstress_lift_scheme:",&
     94                   windstress_lift_scheme
     95        call abort_physic(rname,"bad windstress_lift_scheme value",1)
     96      endif
     97     
     98      firstcall=.false.
     99    endif ! of if (firstcall)
     100
     101
     102    if (windstress_lift_scheme==0) then
     103
     104!   ---------------------------------
     105!   1. Compute horizontal flux: fhoriz
     106!   ---------------------------------
    82107
    83108      do ig=1,ngrid
    84           fhoriz(ig) = 0.      ! initialisation
     109        fhoriz(ig) = 0.      ! initialisation
    85110
    86 c         Selection of points where surface dust is available
    87 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    88 c         if (latid(ig).ge.80.) goto 99  ! N permanent  polar caps
    89 c         if (latid(ig).le.-80.) goto 99 ! S polar deposits
    90 c         if  ((longd(ig).ge.-141. .and. longd(ig).le.-127.)
    91 c    &   .and.(latid(ig).ge.12.   .and. latid(ig).le.23.))goto 99 ! olympus
    92 c         if  ((longd(ig).ge.-125. .and. longd(ig).le.-118.)
    93 c    &   .and.(latid(ig).ge.-12.   .and. latid(ig).le.-6.))goto 99 ! Arsia
    94 c         if  ((longd(ig).ge.-116. .and. longd(ig).le.-109.)
    95 c    &   .and.(latid(ig).ge.-5.   .and. latid(ig).le. 5.))goto 99 ! pavonis
    96 c         if  ((longd(ig).ge.-109. .and. longd(ig).le.-100.)
    97 c    &   .and.(latid(ig).ge. 7.   .and. latid(ig).le. 16.))goto 99 ! ascraeus
    98 c         if  ((longd(ig).ge.  61. .and. longd(ig).le.  63.)
    99 c    &   .and.(latid(ig).ge. 63. .and. latid(ig).le. 64.))goto 99 !weird point
    100           if (co2ice(ig).gt.0.) goto 99
     111!         Selection of points where surface dust is available
     112!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     113!         if (latid(ig).ge.80.) goto 99  ! N permanent  polar caps
     114!         if (latid(ig).le.-80.) goto 99 ! S polar deposits
     115!         if  ((longd(ig).ge.-141. .and. longd(ig).le.-127.)
     116!    &   .and.(latid(ig).ge.12.   .and. latid(ig).le.23.))goto 99 ! olympus
     117!         if  ((longd(ig).ge.-125. .and. longd(ig).le.-118.)
     118!    &   .and.(latid(ig).ge.-12.   .and. latid(ig).le.-6.))goto 99 ! Arsia
     119!         if  ((longd(ig).ge.-116. .and. longd(ig).le.-109.)
     120!    &   .and.(latid(ig).ge.-5.   .and. latid(ig).le. 5.))goto 99 ! pavonis
     121!         if  ((longd(ig).ge.-109. .and. longd(ig).le.-100.)
     122!    &   .and.(latid(ig).ge. 7.   .and. latid(ig).le. 16.))goto 99 ! ascraeus
     123!         if  ((longd(ig).ge.  61. .and. longd(ig).le.  63.)
     124!    &   .and.(latid(ig).ge. 63. .and. latid(ig).le. 64.))goto 99 !weird point
     125        if (co2ice(ig).gt.0.) goto 99
    101126
    102127
    103 c         Is the wind strong enough ?
    104 c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~
    105           ust = sqrt(stress_seuil/rho(ig))
    106           us = pcdh(ig) /  sqrt(pcdh_true(ig)) ! ustar=cd*v /sqrt(cd)
    107           if (us.gt.ust) then
    108 c            If lifting ?
    109 c            Calcul du flux suivant Marticorena ( en fait white (1979))
     128!       Is the wind strong enough ?
     129!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     130        ust = sqrt(stress_seuil/rho(ig))
     131        us = pcdh(ig) /  sqrt(pcdh_true(ig)) ! ustar=cd*v /sqrt(cd)
     132        if (us.gt.ust) then
     133!         If lifting ?
     134!         Compute flux according to Marticorena (in fact white (1979))
    110135
    111              fhoriz(ig) = 2.61*(rho(ig)/g) *
    112      &      (us -ust) * (us + ust)**2
    113           end if
    114  99      continue
    115       end do
     136          fhoriz(ig) = 2.61*(rho(ig)/g) * (us -ust) * (us + ust)**2
     137         endif
    116138
    117 c     -------------------------------------
    118 c     Computing vertical flux and diffusion
    119 c     -------------------------------------
     139 99   continue
     140      enddo ! of do ig=1,ngrid
     141
     142!     -------------------------------------
     143!     2. Compute vertical flux and diffusion
     144!     -------------------------------------
    120145 
    121        do iq=1,nq
    122          do ig=1,ngrid
    123              dqslift(ig,iq)= -alpha_lift(iq)* fhoriz(ig)
     146      do iq=1,nq
     147        do ig=1,ngrid
     148          dqslift(ig,iq)= -alpha_lift(iq)* fhoriz(ig)
    124149
     150!   the vertical flux replaces the turbulent diffusion term which is set to zero
     151!            zb(ig,1) = 0.
     152!c           If surface deposition by turbulence diffusion (impaction...)
     153!c           if(fhoriz(ig).ne.0) then
     154!c           zb(ig,1) = zcdh(ig)*zb0(ig,1)
     155!c           AMount of Surface deposition !
     156!c           pdqs_dif(ig,iq)=pdqs_dif(ig,iq) +
     157!c    &      zb(ig,1)*zq(ig,1,iq)/ptimestep
     158!c          write(*,*) 'zb(1)  = ' ,  zb(ig,1),zcdh(ig),zb0(ig,1)
     159!c
    125160
    126 cc  le  flux vertical remplace le terme de diffusion turb. qui est mis a zero
    127 c            zb(ig,1) = 0.
    128 cc           If surface deposition by turbulence diffusion (impaction...)
    129 cc           if(fhoriz(ig).ne.0) then
    130 cc           zb(ig,1) = zcdh(ig)*zb0(ig,1)
    131 cc           AMount of Surface deposition !
    132 cc           pdqs_dif(ig,iq)=pdqs_dif(ig,iq) +
    133 cc    &      zb(ig,1)*zq(ig,1,iq)/ptimestep
    134 cc          write(*,*) 'zb(1)  = ' ,  zb(ig,1),zcdh(ig),zb0(ig,1)
    135 cc
     161        enddo
     162      enddo
    136163
    137          enddo
    138        enddo
     164    endif ! of if (windstress_lift_scheme==0)
     165   
     166  END SUBROUTINE dust_windstress_lift
    139167
    140       RETURN
    141       END
    142 
     168END MODULE dust_windstress_lift_mod

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

    r3291 r3292  
    1       SUBROUTINE dustdevil(ngrid,nlay,nq, pplev,pu,pv,pt, ptsurf,pq2,
    2      &                pdqdev,pdqs_dev)
    3 
    4       use tracer_mod, only: alpha_devil
    5       use surfdat_h, only: z0_default
    6       USE comcstfi_h
    7       USE mod_phys_lmdz_para, only: is_master,bcast
    8       IMPLICIT NONE
    9 
    10 c=======================================================================
    11 c
    12 c
    13 c  VERSION SPECIAL TRACEURS :
    14 c  Parameterization of dust devil activities
    15 c  to estimate dust lifting
    16 c  The dust devil activity is estimated based on
    17 c  Renno et al. 1998 (JAS 55, 3244-3252) 
    18 c
    19 c  It is proportional to (1-b)*Fs
    20 c
    21 c  With b= [ps**(rcp+1) - ptop**(rcp+1)] / [(ps-ptop)*(rcp+1)* ps**rcp]
    22 c  with ptop pressure of the top of the boundary layer
    23 c       rcp = R/cp
    24 c
    25 c  and Fs the surface sensible heat flux = Cd*|U|*(T(1) -Tsurf)
    26 c       
    27 c=======================================================================
    28 
    29 c-----------------------------------------------------------------------
    30 c   declarations:
    31 c   -------------
    32 
    33 c   arguments:
    34 c   ----------
    35 
    36       INTEGER ngrid,nlay
    37       REAL pplev(ngrid,nlay+1)
    38       REAL pt(ngrid,nlay)
    39       REAL pu(ngrid,nlay)
    40       REAL pv(ngrid,nlay)
    41       REAL pq2(ngrid,nlay+1)
    42       REAL ptsurf(ngrid)
    43 
    44 c    Traceurs :
    45       integer nq
    46       real pdqdev(ngrid,nlay,nq)
    47       real pdqs_dev(ngrid,nq)
     1MODULE dustdevil_mod
     2
     3IMPLICIT NONE
     4
     5INTEGER,SAVE :: dust_devil_scheme ! dust devil scheme number
     6! 0 (default): old Renno et al. 1998 (JAS 55, 3244-3252) scheme
     7
     8!$OMP THREADPRIVATE(dust_devil_scheme)
     9
     10CONTAINS
     11
     12  SUBROUTINE dustdevil(ngrid,nlay,nq, pplev,pu,pv,pt, ptsurf,pq2, &
     13                       pdqdev,pdqs_dev)
     14
     15    use tracer_mod, only: alpha_devil
     16    use surfdat_h, only: z0_default
     17    USE comcstfi_h, ONLY: g, cpp, r, rcp
     18    USE mod_phys_lmdz_para, ONLY: is_master, bcast
     19    USE ioipsl_getin_p_mod, ONLY : getin_p
     20    IMPLICIT NONE
     21
     22!=======================================================================
     23!
     24!  Parameterization of dust devil activities
     25!  to estimate dust lifting
     26!  The dust devil activity is estimated based on
     27!  Renno et al. 1998 (JAS 55, 3244-3252) 
     28!
     29!  It is proportional to (1-b)*Fs
     30!
     31!  With b= [ps**(rcp+1) - ptop**(rcp+1)] / [(ps-ptop)*(rcp+1)* ps**rcp]
     32!  with ptop pressure of the top of the boundary layer
     33!       rcp = R/cp
     34!
     35!  and Fs the surface sensible heat flux = Cd*|U|*(T(1) -Tsurf)
     36!       
     37!=======================================================================
     38
     39!   arguments:
     40!   ----------
     41
     42    INTEGER,INTENT(IN) :: ngrid,nlay
     43    REAL,INTENT(IN) :: pplev(ngrid,nlay+1)
     44    REAL,INTENT(IN) :: pt(ngrid,nlay)
     45    REAL,INTENT(IN) :: pu(ngrid,nlay)
     46    REAL,INTENT(IN) :: pv(ngrid,nlay)
     47    REAL,INTENT(IN) :: pq2(ngrid,nlay+1)
     48    REAL,INTENT(IN) :: ptsurf(ngrid)
     49
     50!    Traceurs :
     51    INTEGER,INTENT(IN) :: nq
     52    real,intent(out) :: pdqdev(ngrid,nlay,nq)
     53    real,intent(out) :: pdqs_dev(ngrid,nq)
    4854     
    49 c   local:
    50 c   ------
    51 
    52       INTEGER ig,l,iq
    53       real Cd, z1
    54       save Cd
    55      
     55!   local:
     56!   ------
     57
     58    INTEGER :: ig,l,iq
     59    real,save :: Cd
    5660!$OMP THREADPRIVATE(Cd)
    57 
    58       LOGICAL firstcall
    59       SAVE firstcall
    60 
     61    real :: z1
     62
     63    LOGICAL,SAVE :: firstcall=.true.
    6164!$OMP THREADPRIVATE(firstcall)
    6265
    63       REAL devila(ngrid)
    64       integer ltop(ngrid)
    65       real b,rho,Fs,wind
    66 
    67 
    68 
    69       REAL  q2top , seuil
    70       SAVE  q2top , seuil
    71       DATA q2top/.5/ ! value of q2 at the top of PBL
    72       DATA seuil/.3/ ! value of minimum dust devil activity for dust lifting
    73      
     66    character(len=20),parameter :: rname="dustdevil"
     67
     68    REAL :: devila(ngrid)
     69    integer :: ltop(ngrid)
     70    real :: b,rho,Fs,wind
     71
     72    REAL,SAVE :: q2top=.5 ! value of q2 at the top of PBL
     73    REAL,SAVE :: seuil=.3 ! value of minimum dust devil activity
     74                          ! for dust lifting     
    7475!$OMP THREADPRIVATE(q2top,seuil)
    7576
    76 
    77       DATA firstcall/.true./
    78 
    79 c   TEMPORAIRE AVEC ANLDEVIL : *************
    80 c        real b_diag(ngrid)
    81 c       real localtime(ngrid)
    82 c       common/temporaire/localtime
    83 c      real ztop(ngrid),magwind(ngrid),t1(ngrid)
    84 c      real rcp ,cpp
    85 c      rcp = kappa
    86 c      cpp = r/rcp
    87 c   **********************************
    88        
    89 
    90 c-----------------------------------------------------------------------
    91 c    initialisation
    92 c    --------------
    93 
    94       ! AS: OK firstcall absolute
     77!-----------------------------------------------------------------------
     78!    initialisation
     79!    --------------
     80
     81    ! AS: OK firstcall absolute
     82
     83    IF (firstcall) THEN
     84
     85      ! get scheme number (default=0, old scheme)
     86      dust_devil_scheme=0
     87      call getin_p("dust_devil_scheme",dust_devil_scheme)
     88
     89      ! sanity check on available dust_devil_scheme values
     90      if (dust_devil_scheme.ne.0) then
     91        write(*,*) trim(rname)//" wrong value for dust_devil_scheme:",&
     92                   dust_devil_scheme
     93        call abort_physic(rname,"bad dust_devil_scheme value",1)
     94      endif
    9595
    9696      if(is_master) then
    97       IF (firstcall) THEN
    9897
    9998        write(*,*) 'In dustdevil :'
    10099        write(*,*) '    q2top= ',q2top,'     seuil= ', seuil
    101100
    102 c A rough estimation of the horizontal drag coefficient Cd
    103 c (the scale heigh is taken to be 13 km since we are typically
    104 c dealing with daytime temperature around 250K.
    105 c
    106          z1= -0.5*13.e3*log(pplev(1,2)/pplev(1,1))
    107          Cd = (0.4/log(z1/z0_default))**2
    108 
    109          firstcall=.false.
    110 
    111 c        Temporaire
    112 c        open(77,file='devil')
     101! A rough estimation of the horizontal drag coefficient Cd
     102! (the scale heigh is taken to be 13 km since we are typically
     103! dealing with daytime temperature around 250K.
     104!
     105        z1= -0.5*13.e3*log(pplev(1,2)/pplev(1,1))
     106        Cd = (0.4/log(z1/z0_default))**2
     107
     108!       Temporaire
     109!       open(77,file='devil')
    113110     
    114       ENDIF
    115       endif !ismaster
    116 
     111      endif !is_master
     112
     113      ! share the info with all cores
    117114      CALL bcast(z1)
    118115      CALL bcast(Cd)
    119       CALL bcast(firstcall)
    120 
    121 c-----------------------------------------------------------------------
    122 c Initialisation
    123       do iq=1,nq
    124        do l=1,nlay
    125            do ig=1,ngrid
    126              pdqdev(ig,l,iq)= 0
    127            end do
    128        end do
     116
     117      firstcall=.false.
     118
     119    ENDIF ! firstcall
     120
     121
     122!-----------------------------------------------------------------------
     123! Initialisation of outputs
     124    do iq=1,nq
     125      do l=1,nlay
     126        do ig=1,ngrid
     127          pdqdev(ig,l,iq)= 0
     128        end do
    129129      end do
    130 
    131 
    132 c-----------------------------------------------------------------------
    133 c      Determining the top of the boundary layer
    134 c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     130    end do
     131    pdqs_dev(1:ngrid,1:nq)=0
     132
     133    if (dust_devil_scheme==0) then
     134!-----------------------------------------------------------------------
     135!      Determining the top of the boundary layer
     136!      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    135137      do ig=1,ngrid
    136          do  l=2,nlay-1
    137             if (pq2(ig,l).lt.q2top)then
    138               ltop(ig)=l
    139               goto 99
    140             end if
    141          enddo
    142  99      continue
    143 
    144 c        ***************************************
    145 cc        if (ptsurf(ig).gt.255)then
    146 c         write(*,*) 'tsurf, ztop (km): ', ig,
    147 c     &   ptsurf(ig), -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
    148 c         if ((ptsurf(ig).gt.50.).and.(
    149 c     &      (-12.*log(pplev(ig,ltop(ig))/pplev(ig,1))).gt.60.))then
    150 c            do l=1,nlay
    151 c             write(*,*) l,
    152 c     &       -12.*log(pplev(ig,l)/pplev(ig,1)),pq2(ig,l)
    153 c            end do
    154 c         end if
    155 cc        end if
    156 c        ***************************************
     138        do  l=2,nlay-1
     139          if (pq2(ig,l).lt.q2top)then
     140            ltop(ig)=l
     141            goto 99
     142          end if
     143        enddo
     144 99     continue
     145
     146!        ***************************************
     147!c        if (ptsurf(ig).gt.255)then
     148!         write(*,*) 'tsurf, ztop (km): ', ig,
     149!     &   ptsurf(ig), -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
     150!         if ((ptsurf(ig).gt.50.).and.(
     151!     &      (-12.*log(pplev(ig,ltop(ig))/pplev(ig,1))).gt.60.))then
     152!            do l=1,nlay
     153!             write(*,*) l,
     154!     &       -12.*log(pplev(ig,l)/pplev(ig,1)),pq2(ig,l)
     155!            end do
     156!         end if
     157!c        end if
     158!        ***************************************
    157159     
    158160      enddo
    159161
    160 c        ***************************************
    161 c        do ig=100,148
    162 c           write(*,*)'time,tsurf,ztop', localtime(ig),ptsurf(ig),
    163 c    &      -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
    164 c        end do
    165 c        ***************************************
    166 
    167 
    168 c   Calculation : dust devil intensity
    169 c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     162!        ***************************************
     163!        do ig=100,148
     164!           write(*,*)'time,tsurf,ztop', localtime(ig),ptsurf(ig),
     165!    &      -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
     166!        end do
     167!        ***************************************
     168
     169
     170!     Calculation : dust devil intensity
     171!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    170172      do ig=1,ngrid
    171173
    172 c --------------------------------------------------
    173 c 1) Version 1 : sensible heat flux using actual wind :
    174 c        Wind magnitude:
    175 c        wind = sqrt(pu(ig,1)**2+pv(ig,1)**2)
    176 c
    177 c --------------------------------------------------
    178 c 2) Version 2 : sensible heat flux using  wind = 15 m/s
    179          wind = 15.
    180 c ----------------------------------------------------
    181 c        Density :
    182          rho=pplev(ig,1)/(R*pt(ig,1))
    183 
    184 c        Sensible heat flux (W.m-2) (>0 if up)
    185          Fs= rho*cpp*Cd * wind
    186      &       * (ptsurf(ig) -pt(ig,1))
    187          b= (pplev(ig,1)**(rcp+1) - pplev(ig,ltop(ig))**(rcp+1)) /
    188      &    ( (pplev(ig,1)-pplev(ig,ltop(ig)))*(rcp+1)*pplev(ig,1)**rcp)
    189 
    190 c        b_diag(ig) = b     ! TEMPORAIRE (pour diagnostique)
    191 
    192 c   Energy flux available to drive dust devil (W.m-2) : (1-b)*Fs
    193 c   Dust devil activity :
    194          devila(ig)= max( 0. , (1-b)*Fs - seuil )
     174! --------------------------------------------------
     175! 1) Version 1 : sensible heat flux using actual wind :
     176!      Wind magnitude:
     177!      wind = sqrt(pu(ig,1)**2+pv(ig,1)**2)
     178!
     179! --------------------------------------------------
     180! 2) Version 2 : sensible heat flux using  wind = 15 m/s
     181        wind = 15.
     182! ----------------------------------------------------
     183!      Density :
     184        rho=pplev(ig,1)/(R*pt(ig,1))
     185
     186!        Sensible heat flux (W.m-2) (>0 if up)
     187        Fs= rho*cpp*Cd * wind * (ptsurf(ig) -pt(ig,1))
     188        b= (pplev(ig,1)**(rcp+1) - pplev(ig,ltop(ig))**(rcp+1)) / &
     189           ( (pplev(ig,1)-pplev(ig,ltop(ig)))*(rcp+1)*pplev(ig,1)**rcp)
     190
     191!     b_diag(ig) = b     ! TEMPORAIRE (pour diagnostique)
     192
     193!   Energy flux available to drive dust devil (W.m-2) : (1-b)*Fs
     194!   Dust devil activity :
     195        devila(ig)= max( 0. , (1-b)*Fs - seuil )
    195196      enddo
    196 c   
    197 c     lifted dust (kg m-2 s-1)  (<0 when lifting)
    198 c     ~~~~~~~~~~ 
     197!   
     198!     lifted dust (kg m-2 s-1)  (<0 when lifting)
     199!     ~~~~~~~~~~ 
    199200      do iq=1,nq
    200          do ig=1,ngrid
    201            pdqs_dev(ig,iq)= - alpha_devil(iq) * devila(ig)
    202          end do
     201        do ig=1,ngrid
     202          pdqs_dev(ig,iq)= - alpha_devil(iq) * devila(ig)
     203        end do
    203204      end do
    204205
    205 c     Injection of dust in the atmosphere (up to the top of pbl)
    206 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     206  !   Injection of dust in the atmosphere (up to the top of pbl)
     207  !   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    207208      do iq=1,nq
    208        do ig=1,ngrid
    209          if (devila(ig).ne.0.) then
    210            do l=1,ltop(ig)
    211              pdqdev(ig,l,iq)=-pdqs_dev(ig,iq)*g/
    212      &        (pplev(ig,1)-pplev(ig,ltop(ig)))
    213            end do
    214          end if
    215        end do
     209        do ig=1,ngrid
     210          if (devila(ig).ne.0.) then
     211            do l=1,ltop(ig)
     212               pdqdev(ig,l,iq)=-pdqs_dev(ig,iq)*g/ &
     213                               (pplev(ig,1)-pplev(ig,ltop(ig)))
     214            end do
     215          end if
     216        end do
    216217      end do
    217 
    218 c *********************************************************     
    219 c     TEMPORAIRE AVEC ANLDEVIL:
    220 c     IF (ngrid.gt.1) THEN
    221 c      do ig=2,ngrid-1
    222 c       write(77,88) lati(ig)*180./pi,localtime(ig),
    223 c    &        -12.*log(pplev(ig,ltop(ig))/pplev(ig,1)),
    224 c    &   devil(ig),min(sqrt(pu(ig,1)**2+pv(ig,1)**2),40.),
    225 c    &   ptsurf(ig)-pt(ig,1),ptsurf(ig),pplev(ig,1)
    226 c      end do   
    227 c88    format (f7.3,1x,f7.3,1x,f6.3,1x,f6.4,1x,f7.4,1x,
    228 c    &        f7.3,1x,f7.3,1x,f9.3)
    229 c      do ig=1,ngrid
    230 c       ztop(ig) = -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
    231 c       magwind(ig) = sqrt(pu(ig,1)**2+pv(ig,1)**2)
    232 c       t1(ig) =max(ptsurf(ig)- pt(ig,1),0.)
    233 c      end do
    234 
    235 c       call write_output('dqs_dev','dqs devil',
    236 c    &               'kg.m-2.s-1',pdqs_dev)
    237 c       call write_output('wind','wind',
    238 c    &               'm.s-1',magwind)
    239 c       call write_output('ztop','top pbl',
    240 c    &               'km',ztop)
    241 c       call write_output('tsurf','tsurf',
    242 c    &               'K',ptsurf)
    243 c       call write_output('T1','T(1)',
    244 c    &               'K',t1)
    245 c       call write_output('b','b',
    246 c    &               ' ',b_diag)
    247 c     END If
    248 c *********************************************************     
    249          
    250       RETURN
    251       END
    252 
    253 
     218   
     219    endif ! of if (dust_devil_scheme==0)
     220     
     221  END SUBROUTINE dustdevil
     222
     223END MODULE dustdevil_mod

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

    r3262 r3292  
    7575     &                          dustscaling_mode, dust_rad_adjust,
    7676     &                          freedust, reff_driven_IRtoVIS_scenario
     77      use dustdevil_mod, only: dustdevil
    7778      use turb_mod, only: q2, wstar, ustar, sensibFlux,
    7879     &                    zmax_th, hfmax_th, turb_resolved

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

    r3285 r3292  
    3838      use vdif_cd_mod, only: vdif_cd
    3939      use lmdz_call_atke, only: call_atke
     40      use dust_windstress_lift_mod, only: dust_windstress_lift
    4041      IMPLICIT NONE
    4142
    … …  
    912913           else
    913914#endif
    914             call dustlift(ngrid,nlay,nq,rho,zcdh_true_tmp,zcdh_tmp,
     915            call dust_windstress_lift(ngrid,nlay,nq,rho,
     916     &                    zcdh_true_tmp,zcdh_tmp,
    915917     &                    pqsurf_tmp(:,igcm_co2),pdqsdif_tmp)
    916918#ifndef MESOSCALE
Note: See TracChangeset for help on using the changeset viewer.