Ignore:
Timestamp:
Jun 18, 2025, 5:12:20 PM (5 weeks ago)
Author:
aborella
Message:

Merge with trunk r5653

Location:
LMDZ6/branches/contrails
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails

  • LMDZ6/branches/contrails/libf/phylmd/Dust/coarsemission.f90

    r5337 r5717  
    55        xlat,xlon,debutphy, &
    66        zu10m,zv10m,wstar,ale_bl,ale_wake, &
     7        nsurfwind,wind10ms,probu, &
    78        scale_param_ssacc,scale_param_sscoa, &
    89        scale_param_dustacc,scale_param_dustcoa, &
    … …  
    5455  REAL, intent(in) ::  xlat(klon)    ! latitudes pour chaque point
    5556  REAL, intent(in) ::  xlon(klon)    ! longitudes pour chaque point
     57  INTEGER, intent(in) ::  nsurfwind
    5658  REAL,DIMENSION(klon),INTENT(IN)    :: zu10m
    5759  REAL,DIMENSION(klon),INTENT(IN)    :: zv10m
    5860  REAL,DIMENSION(klon),INTENT(IN)    :: wstar,Ale_bl,ale_wake
     61  REAL,DIMENSION(klon,nsurfwind),INTENT(IN)    :: wind10ms
     62  REAL,DIMENSION(klon,nsurfwind),INTENT(IN)    :: probu
    5963
    6064  !
    … …  
    190194  param_wstarWAKE(i)=param_wstarWAKEperregion(iregion_wstardust(i))
    191195  ENDDO
    192 
    193 
    194   CALL dustemission( debutphy, xlat, xlon, pctsrf, &
     196 
     197
     198  CALL dustemission( debutphy, xlat, xlon, nsurfwind, pctsrf, &
    195199        zu10m,zv10m,wstar,ale_bl,ale_wake, &
    196200        param_wstarBL, param_wstarWAKE, &
     201        wind10ms, probu, &
    197202        dustsourceacc,dustsourcecoa, &
    198203        dustsourcesco,maskd)

  • LMDZ6/branches/contrails/libf/phylmd/Dust/dustemission_mod.f90

    r5337 r5717  
    1111  INTEGER, PARAMETER     :: nmode=3   ! number of soil-dust modes
    1212  INTEGER, PARAMETER     :: ntyp=5   ! number of soil types
    13   INTEGER, PARAMETER     :: nwb=12   ! number of points for the 10m wind
     13  !INTEGER, PARAMETER     :: nwb=12   ! number of points for the 10m wind
    1414! speed weibull distribution (>=2)
    1515  real   ,parameter     :: z10m=1000. !10m in cm
    … …  
    165165  END SUBROUTINE dustemis_out_init
    166166
    167   SUBROUTINE dustemission( debutphy, xlat, xlon, &    !Input
     167  SUBROUTINE dustemission( debutphy, xlat, xlon, nsurfwind, &    !Input
    168168                          pctsrf,zu10m,zv10m,wstar, & !Input
    169169                          ale_bl,ale_wake, &          !Input
    170                           param_wstarBL, param_wstarWAKE, &  !Input
     170                          param_wstarBL, param_wstarWAKE, & !Input
     171                          wind10ms, probu, & !Input
    171172                          emdustacc,emdustcoa,emdustsco,maskdust)    !Output
    172173  USE dimphy
    … …  
    182183  ! first:
    183184  ! Model grid parameters
     185  INTEGER, INTENT(IN)                :: nsurfwind
    184186  REAL,DIMENSION(klon),     INTENT(IN)     :: xlat
    185187  REAL,DIMENSION(klon),     INTENT(IN)     :: xlon
    … …  
    190192  REAL,DIMENSION(klon),INTENT(IN)          :: ale_bl
    191193  REAL,DIMENSION(klon),INTENT(IN)          :: ale_wake
     194  !REAL,DIMENSION(klon),INTENT(IN)          :: wake_s
     195  !REAL,DIMENSION(klon),INTENT(IN)          :: wake_Cstar
     196  !REAL,DIMENSION(klon),INTENT(IN)          :: zustar
    192197  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarWAKE
    193198  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarBL
    194199 
    195200 
     201  REAL,DIMENSION(klon,nsurfwind), INTENT(IN) :: wind10ms
     202  REAL,DIMENSION(klon,nsurfwind), INTENT(IN) :: probu
     203 
    196204  LOGICAL  :: debutphy ! First physiqs run or not
    197205  ! Intermediate variable: 12 bins emissions
    198   REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: emisbinloc ! vertical emission fluxes
     206  !REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: emisbinloc ! vertical emission fluxes
     207  REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: emisbinloc
    199208
    200209  !OUT variables
    … …  
    206215!  REAL,DIMENSION(klon_glo) :: raux_klon_glo ! auxiliar
    207216
    208 !$OMP THREADPRIVATE(emisbinloc)
     217 INTEGER :: nwb
     218 nwb = nsurfwind
     219!!!$OMP THREADPRIVATE(emisbinloc)
    209220!!!!!!$OMP THREADPRIVATE(maskdust)
    210221  IF (debutphy) THEN
    … …  
    217228
    218229!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
     230  CALL  calcdustemission(debutphy,nsurfwind,zu10m,zv10m,wstar,ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, & !I
     231                         wind10ms,probu,emisbinloc)   !O
    221232
    222233  CALL makemask(maskdust)
    … …  
    654665  varname='A'
    655666  CALL read_surface(varname,Aini)
    656 print *,'beforewritephy',mpi_rank,omp_rank
     667!print *,'beforewritephy',mpi_rank,omp_rank
    657668  CALL writefield_phy("SOLinit",solini,5)
    658669  CALL writefield_phy("Pinit",Pini,5)
    … …  
    662673  CALL writefield_phy("Dinit",Dini,5)
    663674  CALL writefield_phy("Ainit",Aini,5)
    664 print *,'afterwritephy',mpi_rank,omp_rank
     675!print *,'afterwritephy',mpi_rank,omp_rank
    665676
    666677  DO i=1,klon
    … …  
    765776      enddo
    76677730   continue
    767       print*,'IK5'
     778!      print*,'IK5'
    768779      ncl=i-1
    769       print*,'   soil size classes used   ',ncl,' / ',nclass
    770       print*,'   soil size min: ',sizeclass(1),' soil size max: ',sizeclass(ncl)
     780!     print*,'   soil size classes used   ',ncl,' / ',nclass
     781!     print*,'   soil size min: ',sizeclass(1),' soil size max: ',sizeclass(ncl)
    771782      if(ncl.gt.nclass)stop
    772783
    … …  
    775786!if (.true.) then
    776787!c 0: Iversen and White 1982
    777        print *,'Using  Iversen and White 1982 Uth'
     788!        print *,'Using  Iversen and White 1982 Uth'
    778789         do i=1,ncl
    779790            bb=adust*(sizeclass(i)**xdust)+bdust
    … …  
    11071118!--------------------------------------------------------------------------------------
    11081119
    1109   SUBROUTINE calcdustemission(debutphy,zu10m,zv10m,wstar, &
     1120  SUBROUTINE calcdustemission(debutphy,nsurfwind,zu10m,zv10m,wstar, &
    11101121                              ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, &
     1122                              wind10ms, probu, &
    11111123                              emisbin)
    11121124  ! emisions over 12 dust bin
    … …  
    11171129  ! Input
    11181130  LOGICAL, INTENT(IN)                   :: debutphy ! First physiqs run or not
     1131  INTEGER, INTENT(IN)                   :: nsurfwind ! First physiqs run or not
    11191132  REAL,DIMENSION(klon),INTENT(IN)          :: zu10m   ! 10m zonal wind
    11201133  REAL,DIMENSION(klon),INTENT(IN)          :: zv10m   ! meridional 10m wind
    … …  
    11221135  REAL,DIMENSION(klon),INTENT(IN)          :: ale_bl
    11231136  REAL,DIMENSION(klon),INTENT(IN)          :: ale_wake
     1137  REAL,DIMENSION(klon,nsurfwind),INTENT(IN)          :: wind10ms
     1138  REAL,DIMENSION(klon,nsurfwind),INTENT(IN)          :: probu
    11241139 
    11251140  ! Local variables
    … …  
    11301145  REAL,DIMENSION(klon), INTENT(IN) :: param_wstarBL
    11311146  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
     1147  !REAL,DIMENSION(:), ALLOCATABLE,SAVE   :: wind10ms   ! 10m wind distribution in m/s
     1148  !REAL,DIMENSION(:), ALLOCATABLE,SAVE   :: wind10cm   ! 10m wind distribution in cm/s
    11341149  REAL,DIMENSION(klon)                  :: zwstar   
    1135   REAL,DIMENSION(nwb)                :: probu
     1150  !REAL,DIMENSION(nwb)                :: probu
    11361151!  REAL, DIMENSION(nmode) :: fluxN,ftN,adN,fdpN,pN,eN ! in the original code N=1,2,3
    11371152  REAL :: flux1,flux2,flux3,ft1,ft2,ft3
    … …  
    11471162  REAL :: dfec1,dfec2,dfec3,t1,t2,t3,p1,p2,p3,dec,ec
    11481163  ! auxiliar counters
    1149   INTEGER                               :: kwb
     1164  INTEGER                               :: kwb, nwb
    11501165  INTEGER                               :: i,j,k,l,n
    11511166  INTEGER  :: kfin,ideb,ifin,kfin2,istep
    … …  
    11551170  !REAL,DIMENSION(:,:), ALLOCATABLE,SAVE  :: emisbin ! vertical emission fluxes in UNITS for the 12 bins
    11561171  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)
     1172  !$OMP THREADPRIVATE(fluxdust)
     1173!!!$OMP THREADPRIVATE(wind10ms)
     1174!!!$OMP THREADPRIVATE(wind10cm)
     1175
    11601176
    11611177  !----------------------------------------------------
    … …  
    11651181!   ALLOCATE( emisbin(klon,nbins) )
    11661182   ALLOCATE( fluxdust(klon,nmode) )
    1167    ALLOCATE( wind10ms(nwb) )
    1168    ALLOCATE( wind10cm(nwb) )
     1183  ! ALLOCATE( wind10ms(klon,nsurfwind) )
     1184   !ALLOCATE( wind10cm(nwb) )
    11691185  ENDIF !debutphy
    11701186
    … …  
    11901206  !
    11911207    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.
     1208   !  zwstar(i)=sqrt(2.*(param_wstarBL(i)*ale_bl(i)+param_wstarWAKE(i)*ale_wake(i)))
     1209     zwstar(i)=sqrt(2.*(param_wstarBL(i)*ale_bl(i)))
    11951210     ! Wind weibull distribution:
    1196 
     1211           nwb = nsurfwind
     1212!          print*,'GGGGGGGGGGGGGGGGGGGGGGGGG nwb=',nwb
    11971213           DO kwb=1,nwb
    11981214                flux1=0.
    … …  
    12041220! lambda=U10mMOD/gamma(1+1/kref)
    12051221! 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.
    12281222             DO n=1,ntyp
    12291223                   ft1=0.
    … …  
    12681262! Cas ou wsta=0.
    12691263                      cdnms=vkarm/(log(z10m/z0salt))
    1270                       modwm=sqrt((wind10ms(kwb)**2)+(1.2*zwstar(i))**2)
     1264                      modwm=sqrt((wind10ms(i,kwb)**2)+(1.2*zwstar(i))**2)
    12711265                      ustarns=cdnms*modwm*100.
    12721266                    ustarsalt=ustarns
    1273 
     1267!                   print*,'LAAAAAAAAAAAAAAAAAA modwm=',modwm
    12741268
    12751269                   IF(ustarsalt.lt.umin/ceff)GOTO 80
    … …  
    13271321             ENDDO !n=1,ntyp
    1328132270 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)
     1323        fluxdust(i,1)=fluxdust(i,1)+flux1*probu(i,kwb)
     1324        fluxdust(i,2)=fluxdust(i,2)+flux2*probu(i,kwb)
     1325        fluxdust(i,3)=fluxdust(i,3)+flux3*probu(i,kwb)
    13321326   ENDDO !kwb=1,nwb
    13331327      m1dflux(i)=10.*fluxdust(i,1)
    … …  
    14101404         enddo
    14111405         if(kfin.ge.nclass)then
    1412             print*,'$$$$ Tables dimension problem:',kfin,'>',nclass
     1406!           print*,'$$$$ Tables dimension problem:',kfin,'>',nclass
    14131407         endif
    14141408!---------------       

  • LMDZ6/branches/contrails/libf/phylmd/Dust/phytracr_spl_mod.F90

    r5618 r5717  
    804804                      beta_fisrt,beta_v1,                              &  ! I
    805805                      zu10m,zv10m,wstar,ale_bl,ale_wake,               &  ! I
     806                      nsurfwind,wind10ms,probu,                        &  ! I 
    806807                      d_tr_dyn,tr_seri)                                            ! O
    807808!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    847848!  divers:
    848849!  -------
    849 !
     850      INTEGER, intent(in) ::  nsurfwind
     851      REAL,DIMENSION(klon,nsurfwind),INTENT(IN)    :: wind10ms
     852      REAL,DIMENSION(klon,nsurfwind),INTENT(IN)    :: probu
    850853      real,intent(in) :: pdtphys  ! pas d'integration pour la physique (seconde)
    851854      REAL, intent(in):: jD_cur, jH_cur
    … …  
    21532156                        rlat,rlon,debutphy,                                &
    21542157                        zu10m,zv10m,wstar,ale_bl,ale_wake,                 &
     2158                        nsurfwind,wind10ms,probu,                          &
    21552159                        scale_param_ssacc,scale_param_sscoa,               &
    21562160                        scale_param_dustacc,scale_param_dustcoa,           &
Note: See TracChangeset for help on using the changeset viewer.