Skip to content
Snippets Groups Projects
Commit b3d6856e authored by RODIER Quentin's avatar RODIER Quentin
Browse files

P.Tulet 18/07/2024: dust+salt sedimentation: remove protection, and move... 

P.Tulet 18/07/2024: dust+salt sedimentation: remove protection, and move computation in time-spltting for better stability with 3 prognostic moments
parent 86c63271
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/MNH/sedim_dust.f90
+ 37
− 60
View file @ b3d6856e
......@@ -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
src/MNH/sedim_salt.f90
+ 38
− 57
View file @ b3d6856e
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Documentation | Mentions légales | CGU | Accessibilité : non conforme