Index: trunk/LMDZ.GENERIC/libf/phystd/blackl.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/blackl.F (revision 955) +++ trunk/LMDZ.GENERIC/libf/phystd/blackl.F (revision 959) @@ -4,5 +4,5 @@ ! physical constants - sigma=5.6693d-08 + sigma=5.67032D-8 pi=datan(1.d0)*4.d0 c0=2.9979d+08 Index: trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F (revision 955) +++ trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F (revision 959) @@ -75,14 +75,12 @@ ! 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,12 +95,10 @@ 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,6 +107,6 @@ 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,5 +133,5 @@ 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,5 +155,5 @@ 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,6 +177,6 @@ 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,7 +195,7 @@ 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 Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (revision 955) +++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (revision 959) @@ -326,5 +326,5 @@ 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,5 +816,5 @@ 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,5 +1189,5 @@ ! 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,5 +1203,5 @@ ! 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,5 +1397,10 @@ 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,7 +1431,7 @@ 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,6 +1770,6 @@ 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 Index: trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90 (revision 955) +++ trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90 (revision 959) @@ -128,5 +128,5 @@ real*8, parameter :: SCALEP = 1.00D+2 - real*8, parameter :: sigma = 5.67032e-8 + real*8, parameter :: sigma = 5.67032D-8 real*8 Cmk Index: trunk/LMDZ.GENERIC/libf/phystd/setspi.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/setspi.F90 (revision 955) +++ trunk/LMDZ.GENERIC/libf/phystd/setspi.F90 (revision 959) @@ -68,5 +68,5 @@ mm=0 - forceEC=.false. + forceEC=.true. planckcheck=.true. Index: trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F (revision 955) +++ trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F (revision 959) @@ -43,13 +43,13 @@ 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,7 +57,7 @@ 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,5 +102,5 @@ 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,5 +111,5 @@ 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,9 +127,9 @@ 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,13 +139,13 @@ 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,5 +162,5 @@ 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