Changeset 4950 for LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90

Timestamp:

May 22, 2024, 3:16:36 PM (5 weeks ago)

Author:

lebasn

Message:

StratAer?: New model version (microphysic, composition routine, code cleaning, new params...)

File:

• LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90 (modified) (12 diffs)

LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90

@@ -26 +26 @@
   USE aerophys
   USE infotrac_phy
-  USE phys_local_var_mod, ONLY: DENSO4, f_r_wet
-
+  USE phys_local_var_mod, ONLY: DENSO4, DENSO4B, f_r_wet, f_r_wetB
+  USE strataer_local_var_mod, ONLY: flag_new_strat_compo
+  
   IMPLICIT NONE
 
@@ -43 +44 @@
   ! local variables in coagulation routine
   INTEGER                                       :: i,j,k,nb,ilon,ilev
-  REAL, DIMENSION(nbtr_bin)                     :: radius ! aerosol particle radius in each bin [m]
+  REAL, DIMENSION(nbtr_bin)                     :: radiusdry ! dry aerosol particle radius in each bin [m]
+  REAL, DIMENSION(nbtr_bin)                     :: radiuswet ! wet aerosol particle radius in each bin [m]
   REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_t ! Concentration Traceur at time t [U/KgA]
   REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
   REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin)   :: ff   ! Volume fraction of intermediate particles
-  REAL, DIMENSION(nbtr_bin)                     :: V    ! Volume of bins
+  REAL, DIMENSION(nbtr_bin)                     :: Vdry ! Volume dry of bins
+  REAL, DIMENSION(nbtr_bin)                     :: Vwet ! Volume wet of bins
   REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: Vij  ! Volume sum of i and j
   REAL                                          :: eta  ! Dynamic viscosity of air
@@ -82 +85 @@
   include "YOMCST.h"
 
-  DO i=1, nbtr_bin
-   radius(i)=mdw(i)/2.
-   V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
-  ENDDO
-
-  DO j=1, nbtr_bin
-  DO i=1, nbtr_bin
-   Vij(i,j)= V(i)+V(j)
-  ENDDO
-  ENDDO
-
+! ff(i,j,k): Volume fraction of Vi,j that is partitioned to each model bin k
+! just need to be calculated in model initialization because mdw(:) size is fixed
+! no need to recalculate radius, Vdry, Vij, and ff every timestep because it is for  
+! dry aerosols
+  DO i=1, nbtr_bin
+     radiusdry(i)=mdw(i)/2.
+     Vdry(i)=radiusdry(i)**3.  !neglecting factor 4*RPI/3
+     Vwet(i)=0.0
+  ENDDO
+
+  DO j=1, nbtr_bin
+     DO i=1, nbtr_bin
+        Vij(i,j)= Vdry(i)+Vdry(j)
+     ENDDO
+  ENDDO
+  
 !--pre-compute the f(i,j,k) from Jacobson equation 13
   ff=0.0 
@@ -100 +108 @@
     IF (k.EQ.1) THEN
       ff(i,j,k)= 0.0
-    ELSEIF (k.GT.1.AND.V(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.V(k)) THEN
+    ELSEIF (k.GT.1.AND.Vdry(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.Vdry(k)) THEN
       ff(i,j,k)= 1.-ff(i,j,k-1)
     ELSEIF (k.EQ.nbtr_bin) THEN
-      IF (Vij(i,j).GE.v(k)) THEN
+      IF (Vij(i,j).GE.Vdry(k)) THEN
         ff(i,j,k)= 1.
       ELSE
         ff(i,j,k)= 0.0
       ENDIF
-    ELSEIF (k.LE.(nbtr_bin-1).AND.V(k).LE.Vij(i,j).AND.Vij(i,j).LT.V(k+1)) THEN
-      ff(i,j,k)= V(k)/Vij(i,j)*(V(k+1)-Vij(i,j))/(V(k+1)-V(k))
+    ELSEIF (k.LE.(nbtr_bin-1).AND.Vdry(k).LE.Vij(i,j).AND.Vij(i,j).LT.Vdry(k+1)) THEN
+      ff(i,j,k)= Vdry(k)/Vij(i,j)*(Vdry(k+1)-Vij(i,j))/(Vdry(k+1)-Vdry(k))
     ENDIF
   ENDDO
   ENDDO
   ENDDO
-
+! End of just need to be calculated at initialization because mdw(:) size is fixed
+  
   DO ilon=1, klon
   DO ilev=1, klev 
@@ -120 +129 @@
   IF (is_strato(ilon,ilev)) THEN
   !compute actual wet particle radius & volume for every grid box
-  DO i=1, nbtr_bin
-   radius(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
-   V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
-  ENDDO
-
+  IF(flag_new_strat_compo) THEN
+     DO i=1, nbtr_bin
+        radiuswet(i)=f_r_wetB(ilon,ilev,i)*mdw(i)/2.
+        Vwet(i)= radiuswet(i)**3.  !neglecting factor 4*RPI/3
+!!      Vwet(i)= Vdry(i)*(f_r_wetB(ilon,ilev,i)**3)
+     ENDDO
+  ELSE
+     DO i=1, nbtr_bin
+        radiuswet(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
+        Vwet(i)= radiuswet(i)**3.  !neglecting factor 4*RPI/3
+!!      Vwet(i)= Vdry(i)*(f_r_wet(ilon,ilev)**3)
+     ENDDO
+  ENDIF
+  
 !--Calculations for the coagulation kernel---------------------------------------------------------
 
@@ -150 +168 @@
   Di=0.0
   DO i=1, nbtr_bin
-   Kn(i)=mnfrpth/radius(i)
-   Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radius(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
+      Kn(i)=mnfrpth/radiuswet(i)
+      Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radiuswet(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
   ENDDO
 
 !--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
   thvelpar=0.0
-  DO i=1, nbtr_bin
-   m_par(i)=4./3.*RPI*radius(i)**3.*DENSO4(ilon,ilev)*1000.
-   thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
-  ENDDO
+  IF(flag_new_strat_compo) THEN
+     DO i=1, nbtr_bin
+        m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4B(ilon,ilev,i)*1000.
+        thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
+     ENDDO
+  ELSE
+     DO i=1, nbtr_bin
+        m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4(ilon,ilev)*1000.
+        thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
+     ENDDO
+  ENDIF
 
 !--pre-compute the particle mean free path mfppar(i) from equation 22
@@ -171 +196 @@
   delta=0.0
   DO i=1, nbtr_bin
-   delta(i)=((2.*radius(i)+mfppar(i))**3.-(4.*radius(i)**2.+mfppar(i)**2.)**1.5)/ &
-           & (6.*radius(i)*mfppar(i))-2.*radius(i)
-  ENDDO
-
+      delta(i)=((2.*radiuswet(i)+mfppar(i))**3.-(4.*radiuswet(i)**2.+mfppar(i)**2.)**1.5)/ &
+           & (6.*radiuswet(i)*mfppar(i))-2.*radiuswet(i)
+  ENDDO
+
+!   beta(i,j): coagulation kernel (rate coefficient) of 2 colliding particles i,j 
 !--pre-compute the beta(i,j) from equation 17 in Jacobson
   num=0.0
@@ -180 +206 @@
   DO i=1, nbtr_bin
 !
-   num=4.*RPI*(radius(i)+radius(j))*(Di(i)+Di(j))
-   denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
-        & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radius(i)+radius(j)))
-   beta(i,j)=num/denom
+     num=4.*RPI*(radiuswet(i)+radiuswet(j))*(Di(i)+Di(j))
+     denom=(radiuswet(i)+radiuswet(j))/(radiuswet(i)+radiuswet(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
+          & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radiuswet(i)+radiuswet(j)))
+     beta(i,j)=num/denom
 !
 !--compute enhancement factor due to van der Waals forces
    IF (ok_vdw .EQ. 0) THEN      !--no enhancement factor
-     Evdw=1.0
+      Evdw=1.0
    ELSEIF (ok_vdw .EQ. 1) THEN  !--E(0) case
-     AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
-     xvdW = LOG(1.+AvdWi)
-     EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
+      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
+      xvdW = LOG(1.+AvdWi)
+      EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
    ELSEIF (ok_vdw .EQ. 2) THEN  !--E(infinity) case
-     AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
-     xvdW = LOG(1.+AvdWi)
-     EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
+      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
+      xvdW = LOG(1.+AvdWi)
+      EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
    ENDIF
 !
@@ -209 +235 @@
   denom=0.0
   DO j=1, nbtr_bin
-  denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
+     !    fraction of coagulation of k and j that is not giving k
+     denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
   ENDDO
 
@@ -219 +246 @@
     num=0.0
     DO j=1, k
-    numi=0.0
-    DO i=1, k-1
-    numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
+       numi=0.0
+       DO i=1, k-1
+!           
+!           see Jacobson: " In order to conserve volume and volume concentration (which
+!           coagulation physically does) while giving up some accuracy in number concentration" 
+!
+!           Coagulation of i and j giving k
+!           with V(i) and then V(j) because it considers i,j and j,i with the double loop
+!
+!           BUT WHY WET VOLUME V(i) in old STRATAER? tracers are already dry aerosols and coagulation 
+!           kernel beta(i,j) accounts for wet aerosols -> reply below
+!
+!             numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
+            numi=numi+ff(i,j,k)*beta(i,j)*Vdry(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
+       ENDDO
+       num=num+numi
     ENDDO
-    num=num+numi
-    ENDDO
 
 !--calculate new concentration of other bins
-    tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/(1.+pdtcoag*denom)/V(k)
+!      tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) )
+    tr_tp1(ilon,ilev,k)=(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*Vdry(k) )
+!
+!       In constant composition (no dependency on aerosol size because no kelvin effect)
+!       V(l)= (f_r_wet(ilon,ilev)**3)*((mdw(l)/2.)**3) = (f_r_wet(ilon,ilev)**3)*Vdry(i)
+!       so numi and num are proportional (f_r_wet(ilon,ilev)**3)
+!       and so 
+!        tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) )
+!                     =(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num_dry)/( (1.+pdtcoag*denom)*Vdry(k) ) 
+!          with num_dry=...beta(i,j)*Vdry(i)*....
+!       so in old STRATAER (.not.flag_new_strat_compo), it was correct
   ENDIF
 
@@ -234 +282 @@
 !--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
   DO i=1, nbtr_bin
-   tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
+     tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
   ENDDO
 
@@ -240 +288 @@
   ENDDO !--end of loop klev
   ENDDO !--end of loop klon
+! *********************************************
 
 END SUBROUTINE COAGULATE