Changeset 3409 for trunk/LMDZ.PLUTO/libf

Timestamp:

Aug 20, 2024, 12:12:40 PM (4 months ago)

Author:

afalco

Message:

Pluto PCM:
Take into account zeros in aerosol optical properties.
AF

Location:

trunk/LMDZ.PLUTO/libf/phypluto

Files:

• dsolver.F (modified) (4 diffs)
• gfluxv.F (modified) (19 diffs)

trunk/LMDZ.PLUTO/libf/phypluto/dsolver.F

@@ -40 +40 @@
 
       L=2*NL
- 
+
 C     ************MIXED COEFFICENTS**********
 C     THIS VERSION AVOIDS SINGULARITIES ASSOC.
@@ -53 +53 @@
 
 C     EVEN TERMS
- 
+
       DO I=2,LM2,2
         N     = N+1
-        AF(I) = (E1(N)+E3(N))*(GAMA(N+1)-1.)       
+        AF(I) = (E1(N)+E3(N))*(GAMA(N+1)-1.)
         BF(I) = (E2(N)+E4(N))*(GAMA(N+1)-1.)
+        IF (BF(I).eq.0) THEN
+            BF(I) = 1e-16
+        END IF
         CF(I) = 2.0*(1.-GAMA(N+1)**2)
         DF(I) = (GAMA(N+1)-1.) * (CPM1(N+1) - CP(N)) +
      *            (1.-GAMA(N+1))* (CM(N)-CMM1(N+1))
       END DO
- 
+
       N   = 0
       LM1 = L-1
@@ -69 +72 @@
         AF(I) = 2.0*(1.-GAMA(N)**2)
         BF(I) = (E1(N)-E3(N))*(1.+GAMA(N+1))
+        IF (BF(I).eq.0) THEN
+            BF(I) = 1e-16
+        END IF
         CF(I) = (E1(N)+E3(N))*(GAMA(N+1)-1.)
         DF(I) = E3(N)*(CPM1(N+1) - CP(N)) + E1(N)*(CM(N) - CMM1(N+1))
       END DO
- 
+
       AF(L) = E1(NL)-RSF*E3(NL)
       BF(L) = E2(NL)-RSF*E4(NL)
+      IF (BF(L).eq.0) THEN
+        BF(L) = 1e-16
+      END IF
       CF(L) = 0.0
       DF(L) = BSURF-CP(NL)+RSF*CM(NL)
- 
+
       CALL DTRIDGL(L,AF,BF,CF,DF,XK)
- 
+
 C     ***UNMIX THE COEFFICIENTS****
 
@@ -98 +107 @@
 
    28 CONTINUE
- 
+
       RETURN
       END

trunk/LMDZ.PLUTO/libf/phypluto/gfluxv.F

@@ -1 +1 @@
       module gfluxv_mod
-      
+
       implicit none
-      
+
       contains
-      
+
       SUBROUTINE GFLUXV(DTDEL,TDEL,TAUCUMIN,WDEL,CDEL,UBAR0,F0PI,RSF,
      *                  BTOP,BSURF,FMIDP,FMIDM,DIFFV,FLUXUP,FLUXDN)
@@ -12 +12 @@
 C  FOR THE VISIBLE  FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
 C  THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
-C  MEASURED FROM THE TOP OF EACH LAYER.  (DTAU) TOP OF EACH LAYER HAS  
+C  MEASURED FROM THE TOP OF EACH LAYER.  (DTAU) TOP OF EACH LAYER HAS
 C  OPTICAL DEPTH TAU(N).IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
 C  HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
 C  FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
-C THIS SUBROUTINE DIFFERS FROM ITS IR COUNTERPART IN THAT HERE WE SOLVE FOR 
+C THIS SUBROUTINE DIFFERS FROM ITS IR COUNTERPART IN THAT HERE WE SOLVE FOR
 C THE FLUXES DIRECTLY USING THE GENERALIZED NOTATION OF MEADOR AND WEAVOR
 C J.A.S., 37, 630-642, 1980.
-C THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY 
+C THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY
 C VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER
 C
@@ -29 +29 @@
 C WDEL(NLEVEL)  = SINGLE SCATTERING ALBEDO
 C CDEL(NLL)     = ASYMMETRY FACTORS, 0=ISOTROPIC
-C UBARV         = AVERAGE ANGLE, 
+C UBARV         = AVERAGE ANGLE,
 C UBAR0         = SOLAR ZENITH ANGLE
 C F0PI          = INCIDENT SOLAR DIRECT BEAM FLUX
@@ -41 +41 @@
 C added Dec 2002
 C DIFFV         = downward diffuse solar flux at the surface
-C 
+C
 !======================================================================!
 
@@ -76 +76 @@
       NAYER  = L_NLAYRAD
       TAUMAX = L_TAUMAX    !Default is 35.0
-      
+
 !  Delta-Eddington Scaling
 
@@ -91 +91 @@
         FACTOR      = 1.0D0 - WDEL(L)*CDEL(L)**2
         W0(L)       = WDEL(L)*(1.0D0-CDEL(L)**2)/FACTOR
+        IF (W0(L).gt.1) THEN
+            W0(L) = 1
+        END IF
         COSBAR(L)   = CDEL(L)/(1.0D0+CDEL(L))
 
@@ -105 +108 @@
       FACTOR        = 1.0D0 - WDEL(L)*CDEL(L)**2
       W0(L)         = WDEL(L)*(1.0D0-CDEL(L)**2)/FACTOR
+      IF (W0(L).gt.1) THEN
+          W0(L) = 1
+      END IF
       COSBAR(L)     = CDEL(L)/(1.0D0+CDEL(L))
       DTAU(L)       = DTDEL(L)*FACTOR
@@ -123 +129 @@
         ALPHA(L)=SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L) ) )
 
-C       SET OF CONSTANTS DETERMINED BY DOM 
+C       SET OF CONSTANTS DETERMINED BY DOM
 
 !     Quadrature method
@@ -148 +154 @@
 c     but the scaling is Eddington?
 
-        LAMDA(L) = SQRT(G1(L)**2 - G2(L)**2)
+        IF (G1(L) - G2(L) < 1E-15) THEN
+            LAMDA(L) = 0
+        ELSE
+            LAMDA(L) = SQRT(G1(L)**2 - G2(L)**2)
+        END IF
         GAMA(L)  = (G1(L)-LAMDA(L))/G2(L)
       END DO
@@ -156 +166 @@
         G4    = 1.0-G3(L)
         DENOM = LAMDA(L)**2 - 1./UBAR0**2
- 
+
 C       THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
-C       THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN 
+C       THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
 C       THE SCATTERING GOES TO ZERO
 C       PREVENT THIS WITH AN IF STATEMENT
@@ -172 +182 @@
 C       CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
 C       AT THE TOP OF THE LAYER, THAT IS LOWER   OPTICAL DEPTH TAU(L)
- 
+
         CPM1(L) = AP*EXP(-TAU(L)/UBAR0)
         CMM1(L) = AM*EXP(-TAU(L)/UBAR0)
@@ -185 +195 @@
 
 
- 
+
 C     NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
 C     FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
@@ -204 +214 @@
 
 C     NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
- 
+
       DO L=1,L_NLAYRAD-1
         EXPTRM = MIN(TAUMAX,LAMDA(L)*(TAUCUM(2*L+1)-TAUCUM(2*L)))
- 
+
         EP = EXP(EXPTRM)
 
@@ -215 +225 @@
 
 C       THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
-C       THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN 
+C       THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
 C       THE SCATTERING GOES TO ZERO
 C       PREVENT THIS WITH A IF STATEMENT
@@ -237 +247 @@
         FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
         FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
- 
+
 C       ADD THE DIRECT FLUX TO THE DOWNWELLING TERM
 
         FMIDM(L)= FMIDM(L)+UBAR0*F0PI*EXP(-MIN(TAUMID,TAUMAX)/UBAR0)
-   
-      END DO
- 
+
+      END DO
+
 C     FLUX AT THE Ptop layer
 
@@ -256 +266 @@
 
 C     THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
-C     THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN 
+C     THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
 C     THE SCATTERING GOES TO ZERO
 C     PREVENT THIS WITH A IF STATEMENT
@@ -284 +294 @@
 !      fluxdn = fluxdn+UBAR0*F0PI*EXP(-MIN(TAUCUM(1),TAUMAX)/UBAR0)
 !JL18 the line above assumed that the top of the radiative model was P=0
-!   it seems to be better for the IR to use the middle of the last physical layer as the radiative top. 
+!   it seems to be better for the IR to use the middle of the last physical layer as the radiative top.
 !   so we correct the downwelling flux below for the calculation of the heating rate
       fluxdn = fluxdn+UBAR0*F0PI*EXP(-TAUCUM(2)/UBAR0)
@@ -291 +301 @@
 C     DTAU instead of DTAU/2.
 
-      L     = L_NLAYRAD 
+      L     = L_NLAYRAD
       EXPTRM = MIN(TAUMAX,LAMDA(L)*(TAUCUM(L_LEVELS)-
      *                                 TAUCUM(L_LEVELS-1)))
@@ -302 +312 @@
 
 C     THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
-C     THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN 
+C     THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
 C     THE SCATTERING GOES TO ZERO
 C     PREVENT THIS WITH A IF STATEMENT