diff --git a/src/MNH/sedim_dust.f90 b/src/MNH/sedim_dust.f90 index 339f4cfa757042e35c40fd7ec2f3dca87a620469..da27e4869f39a9f9526077c85ce7465e7cab62f5 100644 --- a/src/MNH/sedim_dust.f90 +++ b/src/MNH/sedim_dust.f90 @@ -168,86 +168,67 @@ CALL DUST_VELGRAV(ZSIG(:,:,1:ILU,:), ZRG(:,:,1:ILU,:), & ! ZH=9999. ZVSNUMMAX(:,:,:) = 0. +ZFLUXSED(:,:,:,:) = 0. DO JK=1,ILU ZH(:,:,JK)=PZZ(:,:,JK+1)-PZZ(:,:,JK) - ! Maximum velocity - ZVSNUMMAX(:,:,JK) = MIN(10. *ZH(:,:,JK) / PDTMONITOR,20.) + ! Maximum velocity + ZVSNUMMAX(:,:,JK) = ZH(:,:,JK) / PDTMONITOR ENDDO ! ZHMIN=MINVAL(ZH(:,:,1:ILU)) -ZFLUXSED(:,:,:,:) = 0. -ZFLUXMAX(:,:,:,:) = 0. - !Get loss rate [1/s] instead of [m/s] -DO JN=1,NMODE_DST*3 - !ZVGK(:,:,:,JN)=ZVGK(:,:,:,JN)/ZH(:,:,:) - ZVGK(:,:,1:ILU,JN) = MIN(ZVGK(:,:,1:ILU,JN), ZVSNUMMAX(:,:,1:ILU)) - ZVSMAX=MAXVAL(ZVGK(:,:,1:ILU,JN)) - ZVSMAX=MAX(ZVSMAX,5.) - ! +DO JN=1,NMODE_DST + + ZVSMAX=MAXVAL(ZVGK(:,:,1:ILU,NM3(JN))) ISPLITA = INT(ZVSMAX*PDTMONITOR/ZHMIN)+1 - ISPLITA = MIN(20, ISPLITA) - ! ZTSPLITR = PDTMONITOR / REAL(ISPLITA) - ! - ZFLUXSED(:,:,ILU+1,JN) = 0. - + ZFLUXSED(:,:,ILU+1,NM3(JN)) = 0. + DO JT=1,ISPLITA - ZFLUXSED(:,:,1:ILU,JN)= ZVGK(:,:,1:ILU,JN)* ZPM(:,:,1:ILU,JN) + ZFLUXSED(:,:,1:ILU,NM3(JN))= ZVGK(:,:,1:ILU,NM3(JN))* ZPM(:,:,1:ILU,NM3(JN)) DO JK=1,ILU - ZW(:,:,JK) = ZTSPLITR /(PZZ(:,:,JK+1)-PZZ(:,:,JK)) + ZW(:,:,JK) = ZTSPLITR /(PZZ(:,:,JK+1)-PZZ(:,:,JK)) + ZPM(:,:,JK,NM3(JN))= ZPM(:,:,JK,NM3(JN)) + & + ZW(:,:,JK)*(ZFLUXSED(:,:,JK+1,NM3(JN))- ZFLUXSED(:,:,JK,NM3(JN))) END DO - ! - ! -! - DO JK=1,ILU - ZPM(:,:,JK,JN)= ZPM(:,:,JK,JN) + & - ZW(:,:,JK)*(ZFLUXSED(:,:,JK+1,JN)- ZFLUXSED(:,:,JK,JN)) + ENDDO +ENDDO + +IF (.NOT.(LRGFIX_DST)) THEN ! 2 or 3 moments +DO JN=1,NMODE_DST + ZVSMAX=MAXVAL(ZVGK(:,:,1:ILU,NM0(JN))) + ISPLITA = INT(ZVSMAX*PDTMONITOR/ZHMIN)+1 + ZTSPLITR = PDTMONITOR / REAL(ISPLITA) + ZFLUXSED(:,:,ILU+1,NM0(JN)) = 0. + DO JT=1,ISPLITA + ZFLUXSED(:,:,1:ILU,NM0(JN))= ZVGK(:,:,1:ILU,NM0(JN))* ZPM(:,:,1:ILU,NM0(JN)) + DO JK=1,ILU + ZW(:,:,JK) = ZTSPLITR /(PZZ(:,:,JK+1)-PZZ(:,:,JK)) + ZPM(:,:,JK,NM0(JN))= ZPM(:,:,JK,NM0(JN)) + & + ZW(:,:,JK)*(ZFLUXSED(:,:,JK+1,NM0(JN))- ZFLUXSED(:,:,JK,NM0(JN))) END DO ENDDO - ENDDO -! -DO JN=1,NMODE_DST -!! Calcul pour maintenir le rayon fixe pour la sedimentation : - IF (LRGFIX_DST) THEN - ZPM(:,:,:,NM3(JN)) = ZPM(:,:,:,NM0(JN)) *& - (ZRG(:,:,:,JN)**3)*EXP(4.5 * LOG(ZSIG(:,:,:,JN))**2) - ENDIF -! Calcul pour maintenir la dispersion fixe pour la sedimentation : -! sinon Rg augmente lors de la reconstruction d'une loi log-normale -! a partir des nouveau Mk (sigma diminue plus vite que Rg) -! calcul de M6 en conservant sigma +END IF + +IF (LVARSIG) THEN ! 3 moments ZPM(:,:,:, NM6(JN)) = ZPM(:,:,:,NM0(JN)) & * ( (ZPM(:,:,:,NM3(JN))/ZPM(:,:,:,NM0(JN)))**(1./3.) & * exp(-(3./2.)*LOG(ZSIG(:,:,:,JN))**2))**6 & * exp(18.*LOG(ZSIG(:,:,:,JN))**2) - - IF ((LSEDFIX).AND.(LVARSIG)) THEN - ! calcul de M6 en conservant Rg - - ZPM(:,:,:,NM6(JN)) = ZPM(:,:,:,NM3(JN)) ** 4 / & - (ZRG(:,:,:,JN)**6 * ZPM(:,:,:,NM0(JN))**3) - - END IF -ENDDO -! -IF (LSEDFIX) THEN - LSEDFIX = .FALSE. -ELSE - LSEDFIX = .TRUE. END IF + ! !* 5. Return to concentration in ppv (#/molec_{air}) ! DO JN=1,NMODE_DST IF (LVARSIG) THEN -WHERE ((ZPM(:,:,:,NM0(JN)) .GT. 100.*ZPMIN(NM0(JN))).AND.& - (ZPM(:,:,:,NM3(JN)) .GT. 100.*ZPMIN(NM3(JN))).AND.& - (ZPM(:,:,:,NM6(JN)) .GT. 100.*ZPMIN(NM6(JN)))) +WHERE ((ZPM(:,:,:,NM0(JN)) .GT. 10.*ZPMIN(NM0(JN))).AND.& + (ZPM(:,:,:,NM3(JN)) .GT. 10.*ZPMIN(NM3(JN))).AND.& + (ZPM(:,:,:,NM6(JN)) .GT. 10.*ZPMIN(NM6(JN)))) PSVT(:,:,:,1+(JN-1)*3) = ZPM(:,:,:,NM0(JN)) * XMD / & (XAVOGADRO*PRHODREF(:,:,:)) PSVT(:,:,:,2+(JN-1)*3) = ZPM(:,:,:,NM3(JN)) * XMD*XPI * 4./3.*ZRHOI / & @@ -264,7 +245,7 @@ ELSEWHERE ENDWHERE ELSE IF (LRGFIX_DST) THEN -WHERE ((ZPM(:,:,:,NM3(JN)) .GT. 100*ZPMIN(NM3(JN)))) +WHERE ((ZPM(:,:,:,NM3(JN)) .GT. 10*ZPMIN(NM3(JN)))) PSVT(:,:,:,JN) = ZPM(:,:,:,NM3(JN)) * XMD*XPI * 4./3.*ZRHOI / & (ZMI*PRHODREF(:,:,:)*XM3TOUM3) ELSEWHERE @@ -274,9 +255,8 @@ ENDWHERE ELSE - -WHERE ((ZPM(:,:,:,NM0(JN)) .GT. 100*ZPMIN(NM0(JN))).AND.& - (ZPM(:,:,:,NM3(JN)) .GT. 100*ZPMIN(NM3(JN)))) +WHERE ((ZPM(:,:,:,NM0(JN)) .GT. 10*ZPMIN(NM0(JN))).AND.& + (ZPM(:,:,:,NM3(JN)) .GT. 10*ZPMIN(NM3(JN)))) PSVT(:,:,:,1+(JN-1)*2) = ZPM(:,:,:,NM0(JN)) * XMD / & (XAVOGADRO*PRHODREF(:,:,:)) PSVT(:,:,:,2+(JN-1)*2) = ZPM(:,:,:,NM3(JN)) * XMD*XPI * 4./3.*ZRHOI / & @@ -290,8 +270,5 @@ ENDWHERE END IF ENDDO - -CALL PPP2DUST(PSVT, PRHODREF, PSIG3D=ZSIG, PRG3D=ZRG, PM3D=ZPM) - ! END SUBROUTINE SEDIM_DUST diff --git a/src/MNH/sedim_salt.f90 b/src/MNH/sedim_salt.f90 index 7027377414104596a76b390405a941fe023bc382..d7c24268dc30e108a180169bc786792f63a92619 100644 --- a/src/MNH/sedim_salt.f90 +++ b/src/MNH/sedim_salt.f90 @@ -95,7 +95,7 @@ REAL, DIMENSION(NMODE_SLT) :: ZINIRADIUS REAL :: ZRGMIN, ZSIGMIN REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),NMODE_SLT) :: ZRG, ZSIG,ZDPG REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),3*NMODE_SLT) :: ZPM, ZPMOLD -REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)+1,3*NMODE_SLT) :: ZFLUXSED, ZFLUXMAX +REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)+1,3*NMODE_SLT) :: ZFLUXSED REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZH,ZMU, ZW, ZVSNUMMAX REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),3*NMODE_SLT) :: ZVGK, ZDPK REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),NMODE_SLT) :: ZVG @@ -168,86 +168,68 @@ CALL SALT_VELGRAV(ZSIG(:,:,1:ILU,:), ZRG(:,:,1:ILU,:), & ! ZH=9999. ZVSNUMMAX(:,:,:) = 0. +ZFLUXSED(:,:,:,:) = 0. DO JK=1,ILU ZH(:,:,JK)=PZZ(:,:,JK+1)-PZZ(:,:,JK) ! Maximum velocity - ZVSNUMMAX(:,:,JK) = MIN(10. *ZH(:,:,JK) / PDTMONITOR,20.) + ZVSNUMMAX(:,:,JK) = ZH(:,:,JK) / PDTMONITOR ENDDO ! ZHMIN=MINVAL(ZH(:,:,1:ILU)) -ZFLUXSED(:,:,:,:) = 0. -ZFLUXMAX(:,:,:,:) = 0. - !Get loss rate [1/s] instead of [m/s] -DO JN=1,NMODE_SLT*3 - !ZVGK(:,:,:,JN)=ZVGK(:,:,:,JN)/ZH(:,:,:) - ZVGK(:,:,1:ILU,JN) = MIN(ZVGK(:,:,1:ILU,JN), ZVSNUMMAX(:,:,1:ILU)) - ZVSMAX=MAXVAL(ZVGK(:,:,1:ILU,JN)) - ZVSMAX=MAX(ZVSMAX,5.) - ! +DO JN=1,NMODE_SLT + + ZVSMAX=MAXVAL(ZVGK(:,:,1:ILU,NM3(JN))) ISPLITA = INT(ZVSMAX*PDTMONITOR/ZHMIN)+1 - ISPLITA = MIN(20, ISPLITA) - ! ZTSPLITR = PDTMONITOR / REAL(ISPLITA) - ! - ZFLUXSED(:,:,ILU+1,JN) = 0. + ZFLUXSED(:,:,ILU+1,NM3(JN)) = 0. DO JT=1,ISPLITA - ZFLUXSED(:,:,1:ILU,JN)= ZVGK(:,:,1:ILU,JN)* ZPM(:,:,1:ILU,JN) + ZFLUXSED(:,:,1:ILU,NM3(JN))= ZVGK(:,:,1:ILU,NM3(JN))* ZPM(:,:,1:ILU,NM3(JN)) DO JK=1,ILU - ZW(:,:,JK) = ZTSPLITR /(PZZ(:,:,JK+1)-PZZ(:,:,JK)) + ZW(:,:,JK) = ZTSPLITR /(PZZ(:,:,JK+1)-PZZ(:,:,JK)) + ZPM(:,:,JK,NM3(JN))= ZPM(:,:,JK,NM3(JN)) + & + ZW(:,:,JK)*(ZFLUXSED(:,:,JK+1,NM3(JN))- ZFLUXSED(:,:,JK,NM3(JN))) END DO - ! - ! -! - DO JK=1,ILU - ZPM(:,:,JK,JN)= ZPM(:,:,JK,JN) + & - ZW(:,:,JK)*(ZFLUXSED(:,:,JK+1,JN)- ZFLUXSED(:,:,JK,JN)) + ENDDO +ENDDO + +IF (.NOT.(LRGFIX_SLT)) THEN ! 2 or 3 moments +DO JN=1,NMODE_SLT + ZVSMAX=MAXVAL(ZVGK(:,:,1:ILU,NM0(JN))) + ISPLITA = INT(ZVSMAX*PDTMONITOR/ZHMIN)+1 + ZTSPLITR = PDTMONITOR / REAL(ISPLITA) + ZFLUXSED(:,:,ILU+1,NM0(JN)) = 0. + DO JT=1,ISPLITA + ZFLUXSED(:,:,1:ILU,NM0(JN))= ZVGK(:,:,1:ILU,NM0(JN))* ZPM(:,:,1:ILU,NM0(JN)) + DO JK=1,ILU + ZW(:,:,JK) = ZTSPLITR /(PZZ(:,:,JK+1)-PZZ(:,:,JK)) + ZPM(:,:,JK,NM0(JN))= ZPM(:,:,JK,NM0(JN)) + & + ZW(:,:,JK)*(ZFLUXSED(:,:,JK+1,NM0(JN))- ZFLUXSED(:,:,JK,NM0(JN))) END DO ENDDO - ENDDO -! -DO JN=1,NMODE_SLT -!! Calcul pour maintenir le rayon fixe pour la sedimentation : - IF (LRGFIX_SLT) THEN - ZPM(:,:,:,NM3(JN)) = ZPM(:,:,:,NM0(JN)) *& - (ZRG(:,:,:,JN)**3)*EXP(4.5 * LOG(ZSIG(:,:,:,JN))**2) - ENDIF -! Calcul pour maintenir la dispersion fixe pour la sedimentation : -! sinon Rg augmente lors de la reconstruction d'une loi log-normale -! a partir des nouveau Mk (sigma diminue plus vite que Rg) -! calcul de M6 en conservant sigma +END IF + +IF (LVARSIG_SLT) THEN ! 3 moments ZPM(:,:,:, NM6(JN)) = ZPM(:,:,:,NM0(JN)) & * ( (ZPM(:,:,:,NM3(JN))/ZPM(:,:,:,NM0(JN)))**(1./3.) & * exp(-(3./2.)*LOG(ZSIG(:,:,:,JN))**2))**6 & * exp(18.*LOG(ZSIG(:,:,:,JN))**2) +END IF - IF ((LSEDFIX).AND.(LVARSIG_SLT)) THEN - ! calcul de M6 en conservant Rg - - ZPM(:,:,:,NM6(JN)) = ZPM(:,:,:,NM3(JN)) ** 4 / & - (ZRG(:,:,:,JN)**6 * ZPM(:,:,:,NM0(JN))**3) - - END IF -ENDDO ! -IF (LSEDFIX) THEN - LSEDFIX = .FALSE. -ELSE - LSEDFIX = .TRUE. -END IF ! !* 5. Return to concentration in ppv (#/molec_{air}) ! DO JN=1,NMODE_SLT IF (LVARSIG_SLT) THEN -WHERE ((ZPM(:,:,:,NM0(JN)) .GT. 100.*ZPMIN(NM0(JN))).AND.& - (ZPM(:,:,:,NM3(JN)) .GT. 100.*ZPMIN(NM3(JN))).AND.& - (ZPM(:,:,:,NM6(JN)) .GT. 100.*ZPMIN(NM6(JN)))) +WHERE ((ZPM(:,:,:,NM0(JN)) .GT. 10.*ZPMIN(NM0(JN))).AND.& + (ZPM(:,:,:,NM3(JN)) .GT. 10.*ZPMIN(NM3(JN))).AND.& + (ZPM(:,:,:,NM6(JN)) .GT. 10.*ZPMIN(NM6(JN)))) PSVT(:,:,:,1+(JN-1)*3) = ZPM(:,:,:,NM0(JN)) * XMD / & (XAVOGADRO*PRHODREF(:,:,:)) PSVT(:,:,:,2+(JN-1)*3) = ZPM(:,:,:,NM3(JN)) * XMD*XPI * 4./3.*ZRHOI / & @@ -263,7 +245,7 @@ ELSEWHERE (XAVOGADRO*PRHODREF(:,:,:)*1.d-6) ENDWHERE ELSE IF (LRGFIX_SLT) THEN -WHERE ((ZPM(:,:,:,NM3(JN)) .GT. 100*ZPMIN(NM3(JN)))) +WHERE ((ZPM(:,:,:,NM3(JN)) .GT. 10.*ZPMIN(NM3(JN)))) PSVT(:,:,:,JN) = ZPM(:,:,:,NM3(JN)) * XMD*XPI * 4./3.*ZRHOI / & (ZMI*PRHODREF(:,:,:)*XM3TOUM3_SALT) ELSEWHERE @@ -273,9 +255,10 @@ ENDWHERE ELSE - -WHERE ((ZPM(:,:,:,NM0(JN)) .GT. 100*ZPMIN(NM0(JN))).AND.& - (ZPM(:,:,:,NM3(JN)) .GT. 100*ZPMIN(NM3(JN)))) +WHERE ((ZPM(:,:,:,NM0(JN)) .GT. 10.*ZPMIN(NM0(JN))).AND.& + (ZPM(:,:,:,NM3(JN)) .GT. 10.*ZPMIN(NM3(JN))).AND.& + (ZRG(:,:,:,JN) .GT. 0.05*ZINIRADIUS(JN)).AND.& + (ZRG(:,:,:,JN) .LT. 2*ZINIRADIUS(JN))) PSVT(:,:,:,1+(JN-1)*2) = ZPM(:,:,:,NM0(JN)) * XMD / & (XAVOGADRO*PRHODREF(:,:,:)) PSVT(:,:,:,2+(JN-1)*2) = ZPM(:,:,:,NM3(JN)) * XMD*XPI * 4./3.*ZRHOI / & @@ -290,7 +273,5 @@ ENDWHERE END IF ENDDO -CALL PPP2SALT(PSVT, PRHODREF, PSIG3D=ZSIG, PRG3D=ZRG, PM3D=ZPM) - ! END SUBROUTINE SEDIM_SALT