Changeset 961 for trunk/LMDZ.GENERIC/libf

Timestamp:

May 16, 2013, 11:07:35 AM (12 years ago)

Author:

jleconte

Message:

15/05/2013 == JL

• correction in radiative scheme to enforce double precision
• corrects calculation of ISR

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• callcorrk.F90 (modified) (7 diffs)
• gfluxi.F (modified) (1 diff)
• inifis.F (modified) (1 diff)
• optci.F90 (modified) (9 diffs)
• physiq.F90 (modified) (1 diff)
• setspi.F90 (modified) (5 diffs)
• sfluxi.F (modified) (3 diffs)
• sfluxv.F (modified) (5 diffs)
• sugas_corrk.F90 (modified) (4 diffs)

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

@@ -113 +113 @@
 
       REAL*8 tauaero(L_LEVELS+1,naerkind)
-      REAL*8 nfluxtopv,nfluxtopi,nfluxtop
+      REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
       real*8 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2)
       real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2)
@@ -125 +125 @@
       INTEGER icount
 
-      real fluxtoplanet
       real szangle
       logical global1d
@@ -291 +290 @@
 
 !     how much light we get
-      fluxtoplanet=0
       do nw=1,L_NSPECTV
          stel(nw)=stellarf(nw)/(dist_star**2)
-         fluxtoplanet=fluxtoplanet + stel(nw)
       end do
 
@@ -657 +654 @@
          if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
 
-            fluxtoplanet=0.
-
             if((ngrid.eq.1).and.(global1d))then
                do nw=1,L_NSPECTV
                   stel_fract(nw)= stel(nw) * 0.25 / acosz
-                  fluxtoplanet=fluxtoplanet + stel_fract(nw)
                                 ! globally averaged = divide by 4
                                 ! but we correct for solar zenith angle
@@ -669 +663 @@
                do nw=1,L_NSPECTV
                   stel_fract(nw)= stel(nw) * fract(ig)
-                  fluxtoplanet=fluxtoplanet + stel_fract(nw)
                end do
             endif
@@ -679 +672 @@
             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
                  acosz,stel_fract,gweight,                         &
-                 nfluxtopv,nfluxoutv_nu,nfluxgndv_nu,              &
+                 nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,              &
                  fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
 
@@ -709 +702 @@
 
 !     Flux incident at the top of the atmosphere
-         fluxtop_dn(ig)=fluxdnv(1)
+         fluxtop_dn(ig)=fluxtopvdn 
 
          fluxtop_lw(ig)  = real(nfluxtopi)

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

@@ -156 +156 @@
         EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL 
         EM    = 1.0D0/EP
-        TERM  = UBARI/(1.-W0(L)*COSBAR(L))
+        TERM  = UBARI/(1.D0-W0(L)*COSBAR(L))
 
 C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE

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

@@ -549 +549 @@
                write(*,*) "cpp=",cpp
            ENDIF
-         else
-           mugaz=8.314*1000./pr
+!         else
+!           mugaz=8.314*1000./pr
          endif
          call su_gases

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

@@ -120 +120 @@
      ! if we have continuum opacities, we need dz
      if(kastprof)then
-        dz(k) = dpr(k)*(1000.0*8.3145/muvar(k))*TMID(K)/(g*PMID(K))
+        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
         U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
      else
@@ -149 +149 @@
 
 ! continuum absorption
-        DCONT = 0.0 
+        DCONT = 0.0d0 
 
         if(continuum.and.(.not.graybody))then
@@ -163 +163 @@
               endif
 
-              dtemp=0.0
+              dtemp=0.0d0
               if(igas.eq.igas_N2)then
 
@@ -180 +180 @@
                  do jgas=1,ngasmx
                     p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
-                    dtempc  = 0.0
+                    dtempc  = 0.0d0
                     if(jgas.eq.igas_N2)then 
                        interm = indi(nw,igas,jgas)
@@ -286 +286 @@
   end do
 
-
+  DTAUKI(L_LEVELS+1,1:L_NSPECTI,1:L_NGAUSS)=0.d0
   !=======================================================================
   !     Now the full treatment for the layers, where besides the opacity
@@ -315 +315 @@
 
         atemp = 0.
-        if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
+        if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
            do iaer=1,naerkind
               atemp = atemp +                                     &
@@ -324 +324 @@
         else
            WBARI(L,nw,ng) = 0.0D0
-           DTAUI(L,NW,NG) = 1.0E-9
+           DTAUI(L,NW,NG) = 1.0D-9
         endif
 
-        if(btemp(L,nw) .GT. 0.0) then
+        if(btemp(L,nw) .GT. 0.0d0) then
            cosbi(L,NW,NG) = atemp/btemp(L,nw)
         else
@@ -344 +344 @@
               DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
 
-              if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
+              if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
 
                  WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
@@ -350 +350 @@
               else
                  WBARI(L,nw,ng) = 0.0D0
-                 DTAUI(L,NW,NG) = 1.0E-9
+                 DTAUI(L,NW,NG) = 1.0D-9
               endif
 

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

@@ -1433 +1433 @@
          print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
          print*, ISR,ASR,OLR,GND,DYN
+
          if(enertest)then
             print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'

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

@@ -127 +127 @@
 
       do M=1,L_NSPECTI
-         WNOI(M)  = 0.5*(BWNI(M+1)+BWNI(M))
+         WNOI(M)  = 0.5D0*(BWNI(M+1)+BWNI(M))
          DWNI(M)  = BWNI(M+1)-BWNI(M)
-         WAVEI(M) = 1.0E+4/WNOI(M)
-         BLAMI(M) = 0.01/BWNI(M)         
-      end do
-      BLAMI(M) = 0.01/BWNI(M)
+         WAVEI(M) = 1.0D+4/WNOI(M)
+         BLAMI(M) = 0.01D0/BWNI(M)         
+      end do
+      BLAMI(M) = 0.01D0/BWNI(M)
 !     note M=L_NSPECTI+1 after loop due to Fortran bizarreness
 
@@ -148 +148 @@
          a = 1.0D-2/BWNI(NW+1)
          b = 1.0D-2/BWNI(NW)
-         bpa = (b+a)/2.0
-         bma = (b-a)/2.0
+         bpa = (b+a)/2.0D0
+         bma = (b-a)/2.0D0
          do nt=NTstar,NTstop
             T   = dble(NT)/NTfac
@@ -167 +167 @@
          print*,'setspi: Force F=sigma*eps*T^4 for all values of T!'
          do nt=NTstar,NTstop
-            plancksum=0.0
+            plancksum=0.0D0
             T=dble(NT)/NTfac
        
@@ -185 +185 @@
       if(planckcheck)then
          ! check energy conservation at lower temperature boundary
-         plancksum=0.0
+         plancksum=0.0D0
          nt=NTstar
          do NW=1,L_NSPECTI
@@ -195 +195 @@
          
          ! check energy conservation at upper temperature boundary
-         plancksum=0.0
+         plancksum=0.0D0
          nt=NTstop
          do NW=1,L_NSPECTI

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

@@ -70 +70 @@
       NTS   = int(TSURF*NTfac)-NTstar+1
       NTT   = int(TTOP *NTfac)-NTstar+1
-!      NTS   = TSURF*10.0D0-499
-!      NTT   = TTOP *10.0D0-499
 
 !JL12 corrects the surface planck function so that its integral is equal to sigma Tsurf^4
@@ -110 +108 @@
 C         OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
  
-          TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
+!          TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
+          TAUTOP = TAUCUMI(2,NW,NG)
           BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
  
@@ -161 +160 @@
        NG     = L_NGAUSS
 
-       TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
+!       TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
+       TAUTOP = TAUCUMI(2,NW,NG)
        BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
 

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

@@ -1 +1 @@
       SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV,
-     *                  UBAR0,STEL,GWEIGHT,NFLUXTOPV,NFLUXOUTV_nu,
-     *                  NFLUXGNDV_nu,
+     *                  UBAR0,STEL,GWEIGHT,NFLUXTOPV,FLUXTOPVDN,
+     *                  NFLUXOUTV_nu,NFLUXGNDV_nu,
      *                  FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf)
 
@@ -18 +18 @@
       real*8 STEL(L_NSPECTV)
       real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD)
-      real*8 NFLUXTOPV, FLUXUP, FLUXDN
+      real*8 NFLUXTOPV, FLUXUP, FLUXDN,FLUXTOPVDN
       real*8 NFLUXOUTV_nu(L_NSPECTV)
       real*8 NFLUXGNDV_nu(L_NSPECTV)
@@ -37 +37 @@
 
       NFLUXTOPV = 0.0
+      FLUXTOPVDN = 0.0
 
       DO NW=1,L_NSPECTV
@@ -100 +101 @@
           NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*GWEIGHT(NG)*
      *                          (1.0-FZEROV(NW))
+          FLUXTOPVDN = FLUXTOPVDN+FLUXDN*GWEIGHT(NG)*
+     *                          (1.0-FZEROV(NW))
           DO L=1,L_NLAYRAD
             FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*
@@ -156 +159 @@
 
         NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*FZERO
+        FLUXTOPVDN = FLUXTOPVDN+FLUXDN*FZERO
         DO L=1,L_NLAYRAD
           FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*FZERO

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

@@ -372 +372 @@
 
          do nw=1,L_NSPECTI
-            fzeroi(nw) = 0.
+            fzeroi(nw) = 0.d0
 !            do nt=1,L_NTREF
 !               do np=1,L_NPREF
@@ -378 +378 @@
 !                     do ng = 1,L_NGAUSS
 !                        if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
-!                           fzeroi(nw)=fzeroi(nw)+1.
+!                           fzeroi(nw)=fzeroi(nw)+1.d0
 !                        endif
 !                     end do
@@ -388 +388 @@
 
          do nw=1,L_NSPECTV
-            fzerov(nw) = 0.
+            fzerov(nw) = 0.d0
 !            do nt=1,L_NTREF
 !               do np=1,L_NPREF
@@ -394 +394 @@
 !                     do ng = 1,L_NGAUSS
 !                        if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
-!                           fzerov(nw)=fzerov(nw)+1.
+!                           fzerov(nw)=fzerov(nw)+1.d0
 !                        endif
 !                     end do