Ignore:
Timestamp:
Feb 6, 2015, 6:53:57 PM (10 years ago)
Author:
jescribano
Message:

Dust emission scheme changes: (1) Included possibility of use previous dust emission scheme (with 2 dust bins). (2) Parameter of Marticorena and Bergametti 1995 set to it's original value (2.61) for testing purposes with pdtphys=15min. (3) Cleaning ustar calculation.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/dustemission_mod.F90

    r2176 r2196  
    2121  real   , parameter :: dmax=0.2
    2222  integer, parameter :: nspe=nmode*3+1
     23  real   ,parameter     :: vkarm=0.41
     24!JE20150202 : updating scheme to chimere13b <<<
     25! original values
     26!  integer, parameter :: div1=3.
     27!  integer, parameter :: div2=3.
     28!  integer, parameter :: div3=3.
     29!  real   , parameter :: e1=3.61/div1
     30!  real   , parameter :: e2=3.52/div2
     31!  real   , parameter :: e3=3.46/div3
     32!  real   , parameter :: factorflux=10.
     33!  real   , parameter :: rop=2.65 ! particle density g/m3
     34!  real   , parameter :: roa=0.001227  ! air density g/m3
     35!  real   , parameter :: pi=3.14159  !!
     36!  real   , parameter :: gravity=981. !! cm!!
     37!  real   , parameter :: cd=1.*roa/gravity
     38! new values
     39!  logical, parameter :: ok_splatunning=.true.
     40! Div=3 from S. Alfaro (Sow et al ACPD 2011)
     41!JE 20150206
    2342  integer, parameter :: div1=3.
    2443  integer, parameter :: div2=3.
    2544  integer, parameter :: div3=3.
    26 
    27   real   ,parameter     :: vkarm=0.41
    2845  real   , parameter :: e1=3.61/div1
    2946  real   , parameter :: e2=3.52/div2
    3047  real   , parameter :: e3=3.46/div3
     48  real   , parameter :: factorflux=1.
    3149  real   , parameter :: rop=2.65 ! particle density g/m3
    3250  real   , parameter :: roa=0.001227  ! air density g/m3
    3351  real   , parameter :: pi=3.14159  !!
    3452  real   , parameter :: gravity=981. !! cm!!
    35   real   , parameter :: cd=1.*roa/gravity
     53! C=2.61 from Marticorena and Bergametti 1995 instead of Gillete and Chen 2001
     54! (recommended C=1.1  in supply-limited dust source area.. )
     55  real   , parameter :: cd=2.61*roa/gravity
     56!JE20150202>>>>
    3657  real,parameter     :: beta=16300.
    3758  real, parameter, dimension(3) :: diam=(/1.5,6.7,14.2/)
    … …  
    6990!!!  INTEGER,DIMENSION(:),ALLOCATABLE,SAVE :: maskdust
    7091  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: feff
     92  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: feffdbg
    7193
    7294  REAL,DIMENSION(:),  ALLOCATABLE,SAVE :: sizeclass
    … …  
    110132!!!!$OMP THREADPRIVATE(maskdust)
    111133!$OMP THREADPRIVATE(feff)
     134!$OMP THREADPRIVATE(feffdbg)
    112135!$OMP THREADPRIVATE(sizeclass)
    113136!$OMP THREADPRIVATE(sizeclass2)
    … …  
    536559    ALLOCATE( solspe(nats,nspe) )
    537560    ALLOCATE( feff(klon,ntyp) )
     561    ALLOCATE( feffdbg(klon,ntyp) )
    538562    ALLOCATE( sizeclass(nclass) )
    539563    ALLOCATE( sizeclass2(nclass) )
    … …  
    768792!              endif
    769793              feff(i,k)=0.
     794              feffdbg(i,k)=0.
    770795 !         print*,'IKKK A ',i,klon,k,ntyp
    771796            else
    … …  
    776801              bb=log(aeff*(xeff/zos(i,k))**0.8)
    777802              cc=1.-aa/bb
     803              feffdbg(i,k)=cc
    778804       !   print*,'IKKK B1 ',i,klon,k,ntyp
    779805! drag partition between zo1 and zo2
    … …  
    789815              endif
    790816              if (feff(i,k).lt.0.)feff(i,k)=0.
     817              if (feffdbg(i,k).lt.0.)feffdbg(i,k)=0.
    791818              if (feff(i,k).gt.1.)feff(i,k)=1.
     819              if (feffdbg(i,k).gt.1.)feffdbg(i,k)=1.
    792820    !      print*,'IKKK E ',i,klon,k,ntyp
    793821            endif
    … …  
    798826!
    799827    CALL writefield_phy("REPART5",feff(1:klon,1:5),5)
     828    CALL writefield_phy("REPART5dbg",feffdbg(1:klon,1:5),5)
    800829endif
    801830
    … …  
    10331062  REAL :: pdfcum,U10mMOD,pdfu,weilambda
    10341063  REAL :: z0salt,ceff,cerod,cpcent
    1035   REAL :: cdnms,ustarns,modwm,utmin
     1064!JE20150202  !REAL :: cdnms,ustarns,modwm,utmin
     1065  REAL :: cdnms,ustarns,modwm
    10361066  REAL :: fdp1,fdp2,ad1,ad2,ad3,flux_diam
    10371067  REAL :: dfec1,dfec2,dfec3,t1,t2,t3,p1,p2,p3,dec,ec
    … …  
    11201150                   ft2=0.
    11211151                   ft3=0.
    1122 
     1152!JE20150129<<<<
     1153
     1154             IF(.FALSE.) THEN
    11231155!                  nat=int(sol(i,n))
    11241156!                    print *,i,n
    … …  
    11361168             !      if(n.eq.1) print*,'nat1=',nat,'sol1=',sol(i,n)
    11371169                   IF(n.eq.1.and.nat.eq.99)GOTO 80
     1170
     1171             ENDIF
     1172             IF(.TRUE.) THEN
     1173                nat=int(sol(i,n))
     1174                if(n == 1 .and. nat >= 14 .or. nat < 1 .or. nat > 19) GOTO 80
     1175             ENDIF
     1176!JE20150129>>>>
     1177
    11381178                      z0salt=z02(i,n)
    11391179                      ceff=feff(i,n)
    … …  
    11491189                      cdnms=vkarm/(log(z10m/z0salt))
    11501190                      modwm=sqrt((wind10ms(kwb)**2)+(1.2*zwstar(i))**2)
    1151 !                      modwm=sqrt((wind10ms(kwb)**2)+(1.2*sqrt(ale_wake(i))
    1152 !                      )**2)
    1153 !                      modwm=sqrt((wind10ms(kwb)**2)+(1.2*sqrt(Ale_bl(i)) )**2)
    1154 !                      modwm=sqrt((wind10ms(kwb)**2))
    11551191                      ustarns=cdnms*modwm*100.
    1156                       utmin=umin/(cdnms*ceff)
    1157                    IF(wind10cm(kwb).ge.utmin)THEN
    1158                       ustarsalt=ustarns+  &
    1159                     (0.3*(wind10cm(kwb)/100.-utmin/100.)**2.)
    1160                    ELSE
    1161                       ustarsalt=ustarns
    1162                    ENDIF
     1192!JE20150202 <<
     1193! Do not have too much sense.. and is not anymore in the chimere14b version.
     1194!
     1195!                      utmin=umin/(cdnms*ceff)
     1196!                   IF(wind10cm(kwb).ge.utmin)THEN
     1197!                      ustarsalt=ustarns+  &
     1198!                    (0.3*(wind10cm(kwb)/100.-utmin/100.)**2.)
     1199!                   ELSE
     1200!                      ustarsalt=ustarns
     1201!                   ENDIF
     1202! ustarsalt should be :
     1203                    ustarsalt=ustarns
     1204!JE20150202 >>
     1205
     1206
    11631207                   IF(ustarsalt.lt.umin/ceff)GOTO 80
    11641208!                      print*,'ustarsalt = ',ustarsalt
    … …  
    12151259             ENDDO !n=1,ntyp
    1216126070 CONTINUE
     1261!factorflux
    12171262        fluxdust(i,1)=fluxdust(i,1)+flux1*probu(kwb)
    12181263        fluxdust(i,2)=fluxdust(i,2)+flux2*probu(kwb)
    12191264        fluxdust(i,3)=fluxdust(i,3)+flux3*probu(kwb)
    12201265   ENDDO !kwb=1,nwb
    1221       m1dflux(i)=10.*fluxdust(i,1)
    1222       m2dflux(i)=10.*fluxdust(i,2)          ! tous en Kg/m2/s
    1223       m3dflux(i)=10.*fluxdust(i,3)
     1266!JE20150202 <<
     1267!      m1dflux(i)=10.*fluxdust(i,1)
     1268!      m2dflux(i)=10.*fluxdust(i,2)          ! tous en Kg/m2/s
     1269!      m3dflux(i)=10.*fluxdust(i,3)
     1270      m1dflux(i)=factorflux*10.*fluxdust(i,1)
     1271      m2dflux(i)=factorflux*10.*fluxdust(i,2)          ! tous en Kg/m2/s
     1272      m3dflux(i)=factorflux*10.*fluxdust(i,3)
     1273!JE20150202 >>
    12241274
    12251275
    … …  
    12341284   DO i=1,ndistb
    12351285    DO j=1,nbins
    1236     emisbin(k,j) = emisbin(k,j)+10*fluxdust(k,i)*massfrac(i,j) 
     1286!JE20150202 <<
     1287!    emisbin(k,j) = emisbin(k,j)+10*fluxdust(k,i)*massfrac(i,j) 
     1288    emisbin(k,j) = emisbin(k,j)+factorflux*fluxdust(k,i)*massfrac(i,j) 
     1289!JE20150202 >>
    12371290    ENDDO !j, nbind
    12381291   ENDDO  !i, nmode
Note: See TracChangeset for help on using the changeset viewer.