Index: trunk/LMDZ.PLUTO/libf/phypluto/dsolver.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/dsolver.F (revision 3408) +++ trunk/LMDZ.PLUTO/libf/phypluto/dsolver.F (revision 3409) @@ -40,5 +40,5 @@ L=2*NL - + C ************MIXED COEFFICENTS********** C THIS VERSION AVOIDS SINGULARITIES ASSOC. @@ -53,14 +53,17 @@ 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,15 +72,21 @@ 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,5 +107,5 @@ 28 CONTINUE - + RETURN END Index: trunk/LMDZ.PLUTO/libf/phypluto/gfluxv.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/gfluxv.F (revision 3408) +++ trunk/LMDZ.PLUTO/libf/phypluto/gfluxv.F (revision 3409) @@ -1,8 +1,8 @@ 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 +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,5 +29,5 @@ 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,5 +41,5 @@ C added Dec 2002 C DIFFV = downward diffuse solar flux at the surface -C +C !======================================================================! @@ -76,5 +76,5 @@ NAYER = L_NLAYRAD TAUMAX = L_TAUMAX !Default is 35.0 - + ! Delta-Eddington Scaling @@ -91,4 +91,7 @@ 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,4 +108,7 @@ 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,5 +129,5 @@ 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,5 +154,9 @@ 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,7 +166,7 @@ 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,5 +182,5 @@ 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,5 +195,5 @@ - + C NOW CALCULATE THE EXPONENTIAL TERMS NEEDED C FOR THE TRIDIAGONAL ROTATED LAYERED METHOD @@ -204,8 +214,8 @@ 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,5 +225,5 @@ 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,11 +247,11 @@ 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,5 +266,5 @@ 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,5 +294,5 @@ ! 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,5 +301,5 @@ 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,5 +312,5 @@ 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