Changeset 5635 for LMDZ6/trunk/libf/phylmd/Dust/dustemission_mod.f90
- Timestamp:
- Apr 30, 2025, 8:10:21 AM (2 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified LMDZ6/trunk/libf/phylmd/Dust/dustemission_mod.f90 ¶
r5337 r5635 1 1 MODULE dustemission_mod 2 2 3 USE lmdz_surf_wind 3 4 IMPLICIT NONE 4 5 !Parameters … … 11 12 INTEGER, PARAMETER :: nmode=3 ! number of soil-dust modes 12 13 INTEGER, PARAMETER :: ntyp=5 ! number of soil types 13 INTEGER, PARAMETER :: nwb=12 ! number of points for the 10m wind14 !INTEGER, PARAMETER :: nwb=12 ! number of points for the 10m wind 14 15 ! speed weibull distribution (>=2) 15 16 real ,parameter :: z10m=1000. !10m in cm … … 165 166 END SUBROUTINE dustemis_out_init 166 167 167 SUBROUTINE dustemission( debutphy, xlat, xlon, & !Input168 SUBROUTINE dustemission( debutphy, xlat, xlon, nsurfwind, & !Input 168 169 pctsrf,zu10m,zv10m,wstar, & !Input 169 170 ale_bl,ale_wake, & !Input 170 param_wstarBL, param_wstarWAKE, & !Input 171 param_wstarBL, param_wstarWAKE, & !Input 172 wind10ms, probu, & !Input 171 173 emdustacc,emdustcoa,emdustsco,maskdust) !Output 172 174 USE dimphy … … 182 184 ! first: 183 185 ! Model grid parameters 186 INTEGER, INTENT(IN) :: nsurfwind 184 187 REAL,DIMENSION(klon), INTENT(IN) :: xlat 185 188 REAL,DIMENSION(klon), INTENT(IN) :: xlon … … 190 193 REAL,DIMENSION(klon),INTENT(IN) :: ale_bl 191 194 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 192 198 REAL,DIMENSION(klon), INTENT(IN) :: param_wstarWAKE 193 199 REAL,DIMENSION(klon), INTENT(IN) :: param_wstarBL 194 200 195 201 202 REAL,DIMENSION(klon,nsurfwind), INTENT(IN) :: wind10ms 203 REAL,DIMENSION(klon,nsurfwind), INTENT(IN) :: probu 204 196 205 LOGICAL :: debutphy ! First physiqs run or not 197 206 ! 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 199 209 200 210 !OUT variables … … 206 216 ! REAL,DIMENSION(klon_glo) :: raux_klon_glo ! auxiliar 207 217 208 !$OMP THREADPRIVATE(emisbinloc) 218 INTEGER :: nwb 219 nwb = nsurfwind 220 !!!$OMP THREADPRIVATE(emisbinloc) 209 221 !!!!!!$OMP THREADPRIVATE(maskdust) 210 222 IF (debutphy) THEN … … 217 229 218 230 !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, & !I220 emisbinloc) !O231 CALL calcdustemission(debutphy,nsurfwind,zu10m,zv10m,wstar,ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, & !I 232 wind10ms,probu,emisbinloc) !O 221 233 222 234 CALL makemask(maskdust) … … 654 666 varname='A' 655 667 CALL read_surface(varname,Aini) 656 print *,'beforewritephy',mpi_rank,omp_rank668 !print *,'beforewritephy',mpi_rank,omp_rank 657 669 CALL writefield_phy("SOLinit",solini,5) 658 670 CALL writefield_phy("Pinit",Pini,5) … … 662 674 CALL writefield_phy("Dinit",Dini,5) 663 675 CALL writefield_phy("Ainit",Aini,5) 664 print *,'afterwritephy',mpi_rank,omp_rank676 !print *,'afterwritephy',mpi_rank,omp_rank 665 677 666 678 DO i=1,klon … … 765 777 enddo 766 778 30 continue 767 print*,'IK5'779 ! print*,'IK5' 768 780 ncl=i-1 769 770 781 ! print*,' soil size classes used ',ncl,' / ',nclass 782 ! print*,' soil size min: ',sizeclass(1),' soil size max: ',sizeclass(ncl) 771 783 if(ncl.gt.nclass)stop 772 784 … … 775 787 !if (.true.) then 776 788 !c 0: Iversen and White 1982 777 print *,'Using Iversen and White 1982 Uth'789 ! print *,'Using Iversen and White 1982 Uth' 778 790 do i=1,ncl 779 791 bb=adust*(sizeclass(i)**xdust)+bdust … … 1107 1119 !-------------------------------------------------------------------------------------- 1108 1120 1109 SUBROUTINE calcdustemission(debutphy, zu10m,zv10m,wstar, &1121 SUBROUTINE calcdustemission(debutphy,nsurfwind,zu10m,zv10m,wstar, & 1110 1122 ale_bl,ale_wake,param_wstarBL,param_wstarWAKE, & 1123 wind10ms, probu, & 1111 1124 emisbin) 1112 1125 ! emisions over 12 dust bin … … 1117 1130 ! Input 1118 1131 LOGICAL, INTENT(IN) :: debutphy ! First physiqs run or not 1132 INTEGER, INTENT(IN) :: nsurfwind ! First physiqs run or not 1119 1133 REAL,DIMENSION(klon),INTENT(IN) :: zu10m ! 10m zonal wind 1120 1134 REAL,DIMENSION(klon),INTENT(IN) :: zv10m ! meridional 10m wind … … 1122 1136 REAL,DIMENSION(klon),INTENT(IN) :: ale_bl 1123 1137 REAL,DIMENSION(klon),INTENT(IN) :: ale_wake 1138 REAL,DIMENSION(klon,nsurfwind),INTENT(IN) :: wind10ms 1139 REAL,DIMENSION(klon,nsurfwind),INTENT(IN) :: probu 1124 1140 1125 1141 ! Local variables … … 1130 1146 REAL,DIMENSION(klon), INTENT(IN) :: param_wstarBL 1131 1147 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/s1133 REAL,DIMENSION(:), ALLOCATABLE,SAVE :: wind10cm ! 10m wind distribution in cm/s1148 !REAL,DIMENSION(:), ALLOCATABLE,SAVE :: wind10ms ! 10m wind distribution in m/s 1149 !REAL,DIMENSION(:), ALLOCATABLE,SAVE :: wind10cm ! 10m wind distribution in cm/s 1134 1150 REAL,DIMENSION(klon) :: zwstar 1135 REAL,DIMENSION(nwb) :: probu1151 !REAL,DIMENSION(nwb) :: probu 1136 1152 ! REAL, DIMENSION(nmode) :: fluxN,ftN,adN,fdpN,pN,eN ! in the original code N=1,2,3 1137 1153 REAL :: flux1,flux2,flux3,ft1,ft2,ft3 … … 1147 1163 REAL :: dfec1,dfec2,dfec3,t1,t2,t3,p1,p2,p3,dec,ec 1148 1164 ! auxiliar counters 1149 INTEGER :: kwb 1165 INTEGER :: kwb, nwb 1150 1166 INTEGER :: i,j,k,l,n 1151 1167 INTEGER :: kfin,ideb,ifin,kfin2,istep … … 1155 1171 !REAL,DIMENSION(:,:), ALLOCATABLE,SAVE :: emisbin ! vertical emission fluxes in UNITS for the 12 bins 1156 1172 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 1160 1177 1161 1178 !---------------------------------------------------- … … 1165 1182 ! ALLOCATE( emisbin(klon,nbins) ) 1166 1183 ALLOCATE( fluxdust(klon,nmode) ) 1167 ALLOCATE( wind10ms(nwb) )1168 ALLOCATE( wind10cm(nwb) )1184 ! ALLOCATE( wind10ms(klon,nsurfwind) ) 1185 !ALLOCATE( wind10cm(nwb) ) 1169 1186 ENDIF !debutphy 1170 1187 … … 1189 1206 !*,'zwstar=sqrt(2.*(',flag_wstarBL,'ale_bl+0.01*(',flag_wstar,'-100)*ale_wake))' 1190 1207 ! 1208 !CALL surf_wind(klon, nsurfwind, zu10m, zv10m, wstar, param_wstarBL, param_wstarWAKE, ale_bl, wind10ms, probu) 1191 1209 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))) 1195 1212 ! Wind weibull distribution: 1196 1213 nwb = nsurfwind 1214 ! print*,'GGGGGGGGGGGGGGGGGGGGGGGGG nwb=',nwb 1197 1215 DO kwb=1,nwb 1198 1216 flux1=0. … … 1204 1222 ! lambda=U10mMOD/gamma(1+1/kref) 1205 1223 ! gamma function estimated with stirling formula 1206 auxreal=1.+1./kref1207 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)THEN1211 wind10ms(kwb)=kwb*2.*U10mMOD/nwb1212 !original1213 ! 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,weilambda1218 !JE20141205>>1219 1220 probu(kwb)=pdfu*2.*U10mMOD/nwb1221 pdfcum=pdfcum+probu(kwb)1222 IF(probu(kwb).le.1.e-2)GOTO 701223 ELSE1224 wind10ms(kwb)=U10mMOD1225 probu(kwb)=1.1226 ENDIF1227 wind10cm(kwb)=wind10ms(kwb)*100.1228 1224 DO n=1,ntyp 1229 1225 ft1=0. … … 1268 1264 ! Cas ou wsta=0. 1269 1265 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) 1271 1267 ustarns=cdnms*modwm*100. 1272 1268 ustarsalt=ustarns 1273 1269 ! print*,'LAAAAAAAAAAAAAAAAAA modwm=',modwm 1274 1270 1275 1271 IF(ustarsalt.lt.umin/ceff)GOTO 80 … … 1327 1323 ENDDO !n=1,ntyp 1328 1324 70 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) 1332 1328 ENDDO !kwb=1,nwb 1333 1329 m1dflux(i)=10.*fluxdust(i,1) … … 1410 1406 enddo 1411 1407 if(kfin.ge.nclass)then 1412 1408 ! print*,'$$$$ Tables dimension problem:',kfin,'>',nclass 1413 1409 endif 1414 1410 !---------------
Note: See TracChangeset
for help on using the changeset viewer.