- Timestamp:
- Apr 5, 2024, 8:34:42 AM (8 months ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 3 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/changelog.txt
r3290 r3292 4594 4594 == 21/03/2024 == Jliu 4595 4595 Found and fixed another small bug in photochemistry related to commit r3289. 4596 4597 == 05/04/2024 == EM 4598 Code cleanup to prepare the addition of new schemes for dust lifting by 4599 wind stress and dust devils. Renamed "dustlift" routine "dust_windstress_lift" 4600 and 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) 1 MODULE dust_windstress_lift_mod 2 3 IMPLICIT NONE 4 5 INTEGER,SAVE :: windstress_lift_scheme ! wind stress lifting scheme number 6 ! 0 (default): old scheme 7 8 !$OMP THREADPRIVATE(windstress_lift_scheme) 9 10 CONTAINS 11 12 SUBROUTINE dust_windstress_lift(ngrid,nlay,nq,rho, & 13 pcdh_true,pcdh,co2ice, & 14 dqslift) 4 15 5 16 #ifndef MESOSCALE 6 17 use tracer_mod, only: alpha_lift, radius 7 18 #else 8 use tracer_mod, only: alpha_lift, radius,9 & igcm_dust_mass, igcm_dust_number,10 &ref_r0,r3n_q19 use tracer_mod, only: alpha_lift, radius, & 20 igcm_dust_mass, igcm_dust_number, & 21 ref_r0,r3n_q 11 22 #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 14 26 15 c=======================================================================16 c 17 cDust lifting by surface winds18 cComputing flux to the middle of the first layer19 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 !======================================================================= 22 34 23 c----------------------------------------------------------------------- 24 c declarations:25 c -------------35 ! 36 ! arguments: 37 ! ---------- 26 38 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) 30 44 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) 37 46 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) 41 54 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" 52 59 53 60 #ifdef MESOSCALE … … 55 62 !!!! AS: ... stress for lifting 56 63 !!!! AS: you have to compile with -DMESOSCALE to do so 57 REALalpha58 REALr0_lift59 INTEGERierr60 REALulim61 OPEN(99,file='stress.def',status='old',form='formatted' 62 .,iostat=ierr)63 64 65 66 67 68 69 70 71 72 73 alpha_lift(igcm_dust_number)=r3n_q*74 &alpha_lift(igcm_dust_mass)/r0_lift**375 76 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 77 84 #endif 78 85 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 ! --------------------------------- 82 107 83 108 do ig=1,ngrid 84 109 fhoriz(ig) = 0. ! initialisation 85 110 86 cSelection of points where surface dust is available87 c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~88 cif (latid(ig).ge.80.) goto 99 ! N permanent polar caps89 cif (latid(ig).le.-80.) goto 99 ! S polar deposits90 cif ((longd(ig).ge.-141. .and. longd(ig).le.-127.)91 c& .and.(latid(ig).ge.12. .and. latid(ig).le.23.))goto 99 ! olympus92 cif ((longd(ig).ge.-125. .and. longd(ig).le.-118.)93 c& .and.(latid(ig).ge.-12. .and. latid(ig).le.-6.))goto 99 ! Arsia94 cif ((longd(ig).ge.-116. .and. longd(ig).le.-109.)95 c& .and.(latid(ig).ge.-5. .and. latid(ig).le. 5.))goto 99 ! pavonis96 cif ((longd(ig).ge.-109. .and. longd(ig).le.-100.)97 c& .and.(latid(ig).ge. 7. .and. latid(ig).le. 16.))goto 99 ! ascraeus98 cif ((longd(ig).ge. 61. .and. longd(ig).le. 63.)99 c& .and.(latid(ig).ge. 63. .and. latid(ig).le. 64.))goto 99 !weird point100 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 101 126 102 127 103 cIs the wind strong enough ?104 c~~~~~~~~~~~~~~~~~~~~~~~~~~~105 106 107 108 cIf 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)) 110 135 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 116 138 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 ! ------------------------------------- 120 145 121 122 123 146 do iq=1,nq 147 do ig=1,ngrid 148 dqslift(ig,iq)= -alpha_lift(iq)* fhoriz(ig) 124 149 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 125 160 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 136 163 137 enddo 138 enddo 164 endif ! of if (windstress_lift_scheme==0) 165 166 END SUBROUTINE dust_windstress_lift 139 167 140 RETURN 141 END 142 168 END 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) 1 MODULE dustdevil_mod 2 3 IMPLICIT NONE 4 5 INTEGER,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 10 CONTAINS 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) 48 54 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 56 60 !$OMP THREADPRIVATE(Cd) 57 58 LOGICAL firstcall 59 SAVE firstcall 60 61 real :: z1 62 63 LOGICAL,SAVE :: firstcall=.true. 61 64 !$OMP THREADPRIVATE(firstcall) 62 65 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 74 75 !$OMP THREADPRIVATE(q2top,seuil) 75 76 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 95 95 96 96 if(is_master) then 97 IF (firstcall) THEN98 97 99 98 write(*,*) 'In dustdevil :' 100 99 write(*,*) ' q2top= ',q2top,' seuil= ', seuil 101 100 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') 113 110 114 ENDIF115 endif !ismaster 116 111 endif !is_master 112 113 ! share the info with all cores 117 114 CALL bcast(z1) 118 115 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 129 129 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 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 135 137 do ig=1,ngrid 136 137 138 139 140 141 142 99 143 144 c***************************************145 cc if (ptsurf(ig).gt.255)then146 cwrite(*,*) 'tsurf, ztop (km): ', ig,147 c& ptsurf(ig), -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))148 cif ((ptsurf(ig).gt.50.).and.(149 c& (-12.*log(pplev(ig,ltop(ig))/pplev(ig,1))).gt.60.))then150 cdo l=1,nlay151 cwrite(*,*) l,152 c& -12.*log(pplev(ig,l)/pplev(ig,1)),pq2(ig,l)153 cend do154 cend if155 cc end if156 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 ! *************************************** 157 159 158 160 enddo 159 161 160 c***************************************161 cdo ig=100,148162 cwrite(*,*)'time,tsurf,ztop', localtime(ig),ptsurf(ig),163 c& -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))164 cend do165 c***************************************166 167 168 cCalculation : dust devil intensity169 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 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 170 172 do ig=1,ngrid 171 173 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 ) 195 196 enddo 196 c197 clifted dust (kg m-2 s-1) (<0 when lifting)198 c~~~~~~~~~~197 ! 198 ! lifted dust (kg m-2 s-1) (<0 when lifting) 199 ! ~~~~~~~~~~ 199 200 do iq=1,nq 200 201 202 201 do ig=1,ngrid 202 pdqs_dev(ig,iq)= - alpha_devil(iq) * devila(ig) 203 end do 203 204 end do 204 205 205 cInjection 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 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 207 208 do iq=1,nq 208 do ig=1,ngrid209 if (devila(ig).ne.0.) then210 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 do214 end if215 end do209 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 216 217 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 223 END MODULE dustdevil_mod -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r3262 r3292 75 75 & dustscaling_mode, dust_rad_adjust, 76 76 & freedust, reff_driven_IRtoVIS_scenario 77 use dustdevil_mod, only: dustdevil 77 78 use turb_mod, only: q2, wstar, ustar, sensibFlux, 78 79 & zmax_th, hfmax_th, turb_resolved -
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
r3285 r3292 38 38 use vdif_cd_mod, only: vdif_cd 39 39 use lmdz_call_atke, only: call_atke 40 use dust_windstress_lift_mod, only: dust_windstress_lift 40 41 IMPLICIT NONE 41 42 … … 912 913 else 913 914 #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, 915 917 & pqsurf_tmp(:,igcm_co2),pdqsdif_tmp) 916 918 #ifndef MESOSCALE
Note: See TracChangeset
for help on using the changeset viewer.