Ignore:
Timestamp:
Apr 30, 2025, 8:10:21 AM (2 months ago)
Author:
fhourdin
Message:

Modified subrid wind for SPLA (Fredho4Lamine)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • TabularUnified LMDZ6/trunk/libf/phylmd/Dust/dustemission_mod.f90

    r5337 r5635  
    11MODULE dustemission_mod
    22
     3  USE lmdz_surf_wind
    34  IMPLICIT NONE
    45  !Parameters   
     
    1112  INTEGER, PARAMETER     :: nmode=3   ! number of soil-dust modes
    1213  INTEGER, PARAMETER     :: ntyp=5   ! number of soil types
    13   INTEGER, PARAMETER     :: nwb=12   ! number of points for the 10m wind
     14  !INTEGER, PARAMETER     :: nwb=12   ! number of points for the 10m wind
    1415! speed weibull distribution (>=2)
    1516  real   ,parameter     :: z10m=1000. !10m in cm
     
    165166  END SUBROUTINE dustemis_out_init
    166167
    167   SUBROUTINE dustemission( debutphy, xlat, xlon, &    !Input
     168  SUBROUTINE dustemission( debutphy, xlat, xlon, nsurfwind, &    !Input
    168169                          pctsrf,zu10m,zv10m,wstar, & !Input
    169170                          ale_bl,ale_wake, &          !Input
    170                           param_wstarBL, param_wstarWAKE, &  !Input
     171                          param_wstarBL, param_wstarWAKE, & !Input
     172                          wind10ms, probu, & !Input
    171173                          emdustacc,emdustcoa,emdustsco,maskdust)    !Output
    172174  USE dimphy
     
    182184  ! first:
    183185  ! Model grid parameters
     186  INTEGER, INTENT(IN)                :: nsurfwind
    184187  REAL,DIMENSION(klon),     INTENT(IN)     :: xlat
    185188  REAL,DIMENSION(klon),     INTENT(IN)     :: xlon
     
    190193  REAL,DIMENSION(klon),INTENT(IN)          :: ale_bl
    191194  REAL,DIMENSION(klon),INTENT(IN)          :: ale_wake
     195  !REAL,DIMENSION(klon),INTENT(IN)          :: wake_s
     196  !REAL,DIMENSION(klon),INTENT(IN)          :: wake_Cstar
     197  !REAL,DIMENSION(klon),INTENT(IN)          :: zustar
    192198  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarWAKE
    193199  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarBL
    194200 
    195201 
     202  REAL,DIMENSION(klon,nsurfwind), INTENT(IN) :: wind10ms
     203  REAL,DIMENSION(klon,nsurfwind), INTENT(IN) :: probu
     204 
    196205  LOGICAL  :: debutphy ! First physiqs run or not
    197206  ! Intermediate variable: 12 bins emissions
    198   REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: emisbinloc ! vertical emission fluxes
     207  !REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: emisbinloc ! vertical emission fluxes
     208  REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: emisbinloc
    199209
    200210  !OUT variables
     
    206216!  REAL,DIMENSION(klon_glo) :: raux_klon_glo ! auxiliar
    207217
    208 !$OMP THREADPRIVATE(emisbinloc)
     218 INTEGER :: nwb
     219 nwb = nsurfwind
     220!!!$OMP THREADPRIVATE(emisbinloc)
    209221!!!!!!$OMP THREADPRIVATE(maskdust)
    210222  IF (debutphy) THEN
     
    217229
    218230!JE20141124  CALL  calcdustemission(debutphy,zu10m,zv10m,wstar,ale_bl,ale_wake,emisbinloc)
    219   CALL  calcdustemission(debutphy,zu10m,zv10m,wstar,ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, & !I
    220                          emisbinloc)   !O
     231  CALL  calcdustemission(debutphy,nsurfwind,zu10m,zv10m,wstar,ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, & !I
     232                         wind10ms,probu,emisbinloc)   !O
    221233
    222234  CALL makemask(maskdust)
     
    654666  varname='A'
    655667  CALL read_surface(varname,Aini)
    656 print *,'beforewritephy',mpi_rank,omp_rank
     668!print *,'beforewritephy',mpi_rank,omp_rank
    657669  CALL writefield_phy("SOLinit",solini,5)
    658670  CALL writefield_phy("Pinit",Pini,5)
     
    662674  CALL writefield_phy("Dinit",Dini,5)
    663675  CALL writefield_phy("Ainit",Aini,5)
    664 print *,'afterwritephy',mpi_rank,omp_rank
     676!print *,'afterwritephy',mpi_rank,omp_rank
    665677
    666678  DO i=1,klon
     
    765777      enddo
    76677830   continue
    767       print*,'IK5'
     779!      print*,'IK5'
    768780      ncl=i-1
    769       print*,'   soil size classes used   ',ncl,' / ',nclass
    770       print*,'   soil size min: ',sizeclass(1),' soil size max: ',sizeclass(ncl)
     781!     print*,'   soil size classes used   ',ncl,' / ',nclass
     782!     print*,'   soil size min: ',sizeclass(1),' soil size max: ',sizeclass(ncl)
    771783      if(ncl.gt.nclass)stop
    772784
     
    775787!if (.true.) then
    776788!c 0: Iversen and White 1982
    777        print *,'Using  Iversen and White 1982 Uth'
     789!        print *,'Using  Iversen and White 1982 Uth'
    778790         do i=1,ncl
    779791            bb=adust*(sizeclass(i)**xdust)+bdust
     
    11071119!--------------------------------------------------------------------------------------
    11081120
    1109   SUBROUTINE calcdustemission(debutphy,zu10m,zv10m,wstar, &
     1121  SUBROUTINE calcdustemission(debutphy,nsurfwind,zu10m,zv10m,wstar, &
    11101122                              ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, &
     1123                              wind10ms, probu, &
    11111124                              emisbin)
    11121125  ! emisions over 12 dust bin
     
    11171130  ! Input
    11181131  LOGICAL, INTENT(IN)                   :: debutphy ! First physiqs run or not
     1132  INTEGER, INTENT(IN)                   :: nsurfwind ! First physiqs run or not
    11191133  REAL,DIMENSION(klon),INTENT(IN)          :: zu10m   ! 10m zonal wind
    11201134  REAL,DIMENSION(klon),INTENT(IN)          :: zv10m   ! meridional 10m wind
     
    11221136  REAL,DIMENSION(klon),INTENT(IN)          :: ale_bl
    11231137  REAL,DIMENSION(klon),INTENT(IN)          :: ale_wake
     1138  REAL,DIMENSION(klon,nsurfwind),INTENT(IN)          :: wind10ms
     1139  REAL,DIMENSION(klon,nsurfwind),INTENT(IN)          :: probu
    11241140 
    11251141  ! Local variables
     
    11301146  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarBL
    11311147  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE :: fluxdust ! horizonal emission fluxes in UNITS for the nmod soil aerosol modes
    1132   REAL,DIMENSION(:), ALLOCATABLE,SAVE   :: wind10ms   ! 10m wind distribution in m/s
    1133   REAL,DIMENSION(:), ALLOCATABLE,SAVE   :: wind10cm   ! 10m wind distribution in cm/s
     1148  !REAL,DIMENSION(:), ALLOCATABLE,SAVE   :: wind10ms   ! 10m wind distribution in m/s
     1149  !REAL,DIMENSION(:), ALLOCATABLE,SAVE   :: wind10cm   ! 10m wind distribution in cm/s
    11341150  REAL,DIMENSION(klon)                  :: zwstar   
    1135   REAL,DIMENSION(nwb)                :: probu
     1151  !REAL,DIMENSION(nwb)                :: probu
    11361152!  REAL, DIMENSION(nmode) :: fluxN,ftN,adN,fdpN,pN,eN ! in the original code N=1,2,3
    11371153  REAL :: flux1,flux2,flux3,ft1,ft2,ft3
     
    11471163  REAL :: dfec1,dfec2,dfec3,t1,t2,t3,p1,p2,p3,dec,ec
    11481164  ! auxiliar counters
    1149   INTEGER                               :: kwb
     1165  INTEGER                               :: kwb, nwb
    11501166  INTEGER                               :: i,j,k,l,n
    11511167  INTEGER  :: kfin,ideb,ifin,kfin2,istep
     
    11551171  !REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: emisbin ! vertical emission fluxes in UNITS for the 12 bins
    11561172  REAL,DIMENSION(klon,nbins)  :: emisbin ! vertical emission fluxes in UNITS for the 12 bins
    1157 !$OMP THREADPRIVATE(fluxdust)
    1158 !$OMP THREADPRIVATE(wind10ms)
    1159 !$OMP THREADPRIVATE(wind10cm)
     1173  !$OMP THREADPRIVATE(fluxdust)
     1174!!!$OMP THREADPRIVATE(wind10ms)
     1175!!!$OMP THREADPRIVATE(wind10cm)
     1176
    11601177
    11611178  !----------------------------------------------------
     
    11651182!   ALLOCATE( emisbin(klon,nbins) )
    11661183   ALLOCATE( fluxdust(klon,nmode) )
    1167    ALLOCATE( wind10ms(nwb) )
    1168    ALLOCATE( wind10cm(nwb) )
     1184  ! ALLOCATE( wind10ms(klon,nsurfwind) )
     1185   !ALLOCATE( wind10cm(nwb) )
    11691186  ENDIF !debutphy
    11701187
     
    11891206!*,'zwstar=sqrt(2.*(',flag_wstarBL,'ale_bl+0.01*(',flag_wstar,'-100)*ale_wake))'
    11901207  !
     1208    !CALL surf_wind(klon, nsurfwind, zu10m, zv10m, wstar, param_wstarBL, param_wstarWAKE, ale_bl, wind10ms, probu)
    11911209    DO i=1,klon  ! main loop
    1192      zwstar(i)=sqrt(2.*(param_wstarBL(i)*ale_bl(i)+param_wstarWAKE(i)*ale_wake(i)))
    1193      U10mMOD=MAX(woff,sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i)))
    1194      pdfcum=0.
     1210   !  zwstar(i)=sqrt(2.*(param_wstarBL(i)*ale_bl(i)+param_wstarWAKE(i)*ale_wake(i)))
     1211     zwstar(i)=sqrt(2.*(param_wstarBL(i)*ale_bl(i)))
    11951212     ! Wind weibull distribution:
    1196 
     1213           nwb = nsurfwind
     1214!          print*,'GGGGGGGGGGGGGGGGGGGGGGGGG nwb=',nwb
    11971215           DO kwb=1,nwb
    11981216                flux1=0.
     
    12041222! lambda=U10mMOD/gamma(1+1/kref)
    12051223! gamma function estimated with stirling formula
    1206                 auxreal=1.+1./kref
    1207                 weilambda = U10mMOD/exp(auxreal*log(auxreal)-auxreal &
    1208                          - 0.5*log(auxreal/(2.*pi))+1./(12.*auxreal) &
    1209                          -1./(360.*(auxreal**3.))+1./(1260.*(auxreal**5.)))
    1210                 IF(nwb.gt.1)THEN
    1211                    wind10ms(kwb)=kwb*2.*U10mMOD/nwb
    1212 !original
    1213 !                   pdfu=(kref/U10mMOD)*(wind10ms(kwb)/U10mMOD)**(kref-1) &
    1214 !                      *exp(-(wind10ms(kwb)/U10mMOD)**kref)
    1215                    pdfu=(kref/weilambda)*(wind10ms(kwb)/weilambda)**(kref-1) &
    1216                       *exp(-(wind10ms(kwb)/weilambda)**kref)
    1217 !                   !print *,'JEdbg  U10mMOD weilambda  ',U10mMOD,weilambda
    1218 !JE20141205>>
    1219 
    1220                    probu(kwb)=pdfu*2.*U10mMOD/nwb
    1221                    pdfcum=pdfcum+probu(kwb)
    1222                       IF(probu(kwb).le.1.e-2)GOTO 70
    1223                 ELSE
    1224                    wind10ms(kwb)=U10mMOD
    1225                    probu(kwb)=1.
    1226                 ENDIF
    1227              wind10cm(kwb)=wind10ms(kwb)*100.
    12281224             DO n=1,ntyp
    12291225                   ft1=0.
     
    12681264! Cas ou wsta=0.
    12691265                      cdnms=vkarm/(log(z10m/z0salt))
    1270                       modwm=sqrt((wind10ms(kwb)**2)+(1.2*zwstar(i))**2)
     1266                      modwm=sqrt((wind10ms(i,kwb)**2)+(1.2*zwstar(i))**2)
    12711267                      ustarns=cdnms*modwm*100.
    12721268                    ustarsalt=ustarns
    1273 
     1269!                   print*,'LAAAAAAAAAAAAAAAAAA modwm=',modwm
    12741270
    12751271                   IF(ustarsalt.lt.umin/ceff)GOTO 80
     
    13271323             ENDDO !n=1,ntyp
    1328132470 CONTINUE
    1329         fluxdust(i,1)=fluxdust(i,1)+flux1*probu(kwb)
    1330         fluxdust(i,2)=fluxdust(i,2)+flux2*probu(kwb)
    1331         fluxdust(i,3)=fluxdust(i,3)+flux3*probu(kwb)
     1325        fluxdust(i,1)=fluxdust(i,1)+flux1*probu(i,kwb)
     1326        fluxdust(i,2)=fluxdust(i,2)+flux2*probu(i,kwb)
     1327        fluxdust(i,3)=fluxdust(i,3)+flux3*probu(i,kwb)
    13321328   ENDDO !kwb=1,nwb
    13331329      m1dflux(i)=10.*fluxdust(i,1)
     
    14101406         enddo
    14111407         if(kfin.ge.nclass)then
    1412             print*,'$$$$ Tables dimension problem:',kfin,'>',nclass
     1408!           print*,'$$$$ Tables dimension problem:',kfin,'>',nclass
    14131409         endif
    14141410!---------------       
Note: See TracChangeset for help on using the changeset viewer.