Changeset 959 for trunk/LMDZ.GENERIC/libf/phystd

Timestamp:

May 15, 2013, 5:39:25 PM (12 years ago)

Author:

jleconte

Message:

15/05/2013 == JL

• correction in radiative scheme to enforce double precision

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• blackl.F (modified) (1 diff)
• gfluxi.F (modified) (7 diffs)
• physiq.F90 (modified) (7 diffs)
• radcommon_h.F90 (modified) (1 diff)
• setspi.F90 (modified) (1 diff)
• sfluxi.F (modified) (7 diffs)

trunk/LMDZ.GENERIC/libf/phystd/blackl.F

@@ -4 +4 @@
 
       ! physical constants
-      sigma=5.6693d-08
+      sigma=5.67032D-8
       pi=datan(1.d0)*4.d0
       c0=2.9979d+08

trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F

@@ -75 +75 @@
 !           endif
 ! Prevent this with an if statement:
-        if (W0(L).eq.1.) then
-           W0(L) = 0.99999
+        if (W0(L).eq.1.D0) then
+           W0(L) = 0.99999D0
         endif
 !----------------------------------------------------
 
-        ALPHA(L) = SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) )
-        LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI
-
-        !NT2   = TLEV(2*L+2)*10.0D0-499
-        !NT    = TLEV(2*L)*10.0D0-499
+        ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
+        LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
+
         NT    = int(TLEV(2*L)*NTfac)   - NTstar+1
         NT2   = int(TLEV(2*L+2)*NTfac) - NTstar+1
@@ -97 +95 @@
 
       if (W0(L).eq.1.) then
-          W0(L) = 0.99999
+          W0(L) = 0.99999D0
       end if
 
-      ALPHA(L) = SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) )
-      LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI
-
-      !NT    = TLEV(2*L+1)*10.0D0-499
-      !NT2   = TLEV(2*L)*10.0D0-499
+      ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
+      LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
+
       NT    = int(TLEV(2*L+1)*NTfac) - NTstar+1
       NT2   = int(TLEV(2*L)*NTfac)   - NTstar+1
@@ -111 +107 @@
 
       DO L=1,L_NLAYRAD
-        GAMA(L) = (1.0-ALPHA(L))/(1.0+ALPHA(L))
-        TERM    = UBARI/(1.0-W0(L)*COSBAR(L))
+        GAMA(L) = (1.0D0-ALPHA(L))/(1.0D0+ALPHA(L))
+        TERM    = UBARI/(1.0D0-W0(L)*COSBAR(L))
 
 C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
@@ -137 +133 @@
 
         EP    = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX))
-        EM    = 1.0/EP
+        EM    = 1.0D0/EP
         E1(L) = EP+GAMA(L)*EM
         E2(L) = EP-GAMA(L)*EM
@@ -159 +155 @@
         DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L)
         EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL 
-        EM    = 1.0/EP
+        EM    = 1.0D0/EP
         TERM  = UBARI/(1.-W0(L)*COSBAR(L))
 
@@ -181 +177 @@
 
       EP   = EXP(MIN((LAMDA(L)*DTAU(L)),TAUMAX)) ! CLIPPED EXPONENTIAL 
-      EM   = 1.0/EP
-      TERM = UBARI/(1.-W0(L)*COSBAR(L))
+      EM   = 1.0D0/EP
+      TERM = UBARI/(1.D0-W0(L)*COSBAR(L))
 
 C     CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
@@ -199 +195 @@
 C     FLUX AT THE PTOP LEVEL
 
-      EP   = 1.0
-      EM   = 1.0
-      TERM = UBARI/(1.0-W0(1)*COSBAR(1))
+      EP   = 1.0D0
+      EM   = 1.0D0
+      TERM = UBARI/(1.0D0-W0(1)*COSBAR(1))
 
 C     CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE

trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

@@ -326 +326 @@
       real mass(ngrid,nlayermx),massarea(ngrid,nlayermx)
       real dEtot, dEtots, AtmToSurf_TurbFlux
-      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
+      real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
       real dEzRadsw(ngrid,nlayermx),dEzRadlw(ngrid,nlayermx),dEzdiff(ngrid,nlayermx)
       real dEdiffs(ngrid),dEdiff(ngrid)
@@ -816 +816 @@
             dEtotSW  = cpp*SUM(massarea(:,:)*zdtsw(:,:))/totarea
             dEtotLW  = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea
-            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea
+            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea !JL13 carefull, albedo can have changed since the last time we called corrk
             dEtotsLW = SUM((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:))/totarea
             dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
@@ -1189 +1189 @@
            ! find qtot
            if(watertest)then
-              iq=3
+              iq=igcm_h2o_ice
               dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
               dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
@@ -1203 +1203 @@
            ! find qtot
            if(watertest)then
-              iq=3
+              iq=igcm_h2o_ice
               dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
               dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
@@ -1397 +1397 @@
            print*,Ts1,Ts2,Ts3,TsS
          endif
-      endif
+      else
+         if (cursor == 1) then
+           print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
+           print*,Ts1,Ts2,Ts3
+         endif         
+      end if
 
 !     ---------------------------------------------------------
@@ -1426 +1431 @@
          endif
 
+         print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
+         print*, ISR,ASR,OLR,GND,DYN
          if(enertest)then
-            print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
-            print*, ISR,ASR,OLR,GND,DYN
             print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
             print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
@@ -1765 +1770 @@
           do ig=1,ngrid
              if(tau_col(ig).gt.1.e3)then
-                print*,'WARNING: tau_col=',tau_col(ig)
-                print*,'at ig=',ig,'in PHYSIQ'
+!                print*,'WARNING: tau_col=',tau_col(ig)
+!                print*,'at ig=',ig,'in PHYSIQ'
              endif
           end do

trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90

@@ -128 +128 @@
       real*8, parameter :: SCALEP = 1.00D+2
 
-      real*8, parameter :: sigma = 5.67032e-8
+      real*8, parameter :: sigma = 5.67032D-8
 
       real*8 Cmk

trunk/LMDZ.GENERIC/libf/phystd/setspi.F90

@@ -68 +68 @@
       mm=0
 
-      forceEC=.false.
+      forceEC=.true.
       planckcheck=.true.
 

trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F

@@ -43 +43 @@
 C     ZERO THE NET FLUXES
     
-      NFLUXTOPI = 0.0
+      NFLUXTOPI = 0.0D0
 
       DO NW=1,L_NSPECTI
-        NFLUXTOPI_nu(NW) = 0.0
+        NFLUXTOPI_nu(NW) = 0.0D0
         DO L=1,L_NLAYRAD
-           FLUXUPI_nu(L,NW) = 0.0
-
-              fup_tmp(nw)=0.0
-              fdn_tmp(nw)=0.0
+           FLUXUPI_nu(L,NW) = 0.0D0
+
+              fup_tmp(nw)=0.0D0
+              fdn_tmp(nw)=0.0D0
 
         END DO
@@ -57 +57 @@
 
       DO L=1,L_NLAYRAD
-        FMNETI(L)  = 0.0
-        FLUXUPI(L) = 0.0
-        FLUXDNI(L) = 0.0
+        FMNETI(L)  = 0.0D0
+        FLUXUPI(L) = 0.0D0
+        FLUXDNI(L) = 0.0D0
       END DO
  
@@ -102 +102 @@
          
           if(TAUGSURF(NW,NG).lt. TLIMIT) then
-            fzero = fzero + (1.0-FZEROI(NW))*GWEIGHT(NG)
+            fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
             goto 30
           end if
@@ -111 +111 @@
  
           TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
-          BTOP   = (1.0-EXP(-TAUTOP/UBARI))*PLTOP
+          BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
  
 C         WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
@@ -127 +127 @@
 
           NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)*
-     *                           (1.0-FZEROI(NW))
+     *                           (1.0D0-FZEROI(NW))
 
 c         and same thing by spectral band... (RDW)
           NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
-     *      +FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW))
+     *      +FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0D0-FZEROI(NW))
 
 
@@ -139 +139 @@
 
             FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*
-     *                              GWEIGHT(NG)*(1.0-FZEROI(NW))
+     *                              GWEIGHT(NG)*(1.0D0-FZEROI(NW))
             FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*
-     *                                (1.0-FZEROI(NW))
+     *                                (1.0D0-FZEROI(NW))
             FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*GWEIGHT(NG)*
-     *                                (1.0-FZEROI(NW))
+     *                                (1.0D0-FZEROI(NW))
 
 c         and same thing by spectral band... (RW)
             FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + 
-     *                FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW))
+     *                FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*(1.0D0-FZEROI(NW))
 
           END DO
@@ -162 +162 @@
 
        TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
-       BTOP   = (1.0-EXP(-TAUTOP/UBARI))*PLTOP
+       BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
 
 C      WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM