Changeset 651


Ignore:
Timestamp:
May 4, 2012, 7:06:09 PM (13 years ago)
Author:
jleconte
Message:
  • Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
  • Minor name changes in watercommon
  • Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

  • Added consistency check in inifis
  • Moved 1d water initialization from physiqu to rcm1d
Location:
trunk/LMDZ.GENERIC
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r650 r651  
    667667- Added consistency check in inifis
    668668- Moved 1d water initialization from physiqu to rcm1d
     669- All enertests in physiq written in a matricial (F90) way. The rest of physiqu should follow soon

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

    r650 r651  
    314314
    315315!     to test energy conservation (RW+JL)
    316       real dEtot, dEtots, masse, AtmToSurf_TurbFlux
     316      real mass(ngridmx,nlayermx),massarea(ngridmx,nlayermx)
     317      real dEtot, dEtots, AtmToSurf_TurbFlux
    317318      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
    318319      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdiff(ngridmx,nlayermx)
    … …  
    626627            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
    627628            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
     629            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/g
     630            massarea(ig,l)=mass(ig,l)*area(ig)
    628631         enddo
    629632      enddo
    … …  
    804807! test energy conservation
    805808         if(enertest)then
    806             dEtotSW  = 0.0
    807             dEtotLW  = 0.0
    808             dEtotsSW = 0.0
    809             dEtotsLW = 0.0
    810             dEzRadsw(:,:)=0.
    811             dEzRadlw(:,:)=0.
    812             do ig = 1, ngrid
    813                do l = 1, nlayer
    814                   masse  = (pplev(ig,l) - pplev(ig,l+1))/g
    815                   dEtotSW = dEtotSW + cpp*masse*zdtsw(ig,l)*area(ig)
    816                   dEtotLW = dEtotLW + cpp*masse*zdtlw(ig,l)*area(ig)
    817                   dEzRadsw(ig,l)=cpp*masse*zdtsw(ig,l)
    818                   dEzRadlw(ig,l)=cpp*masse*zdtlw(ig,l)
    819                enddo
    820                dEtotsSW = dEtotsSW + fluxsurf_sw(ig)*(1.-albedo(ig))*area(ig)
    821                dEtotsLW = dEtotsLW + emis(ig)*fluxsurf_lw(ig)*area(ig)-zplanck(ig)*area(ig)
    822             enddo
    823             dEtotSW  = dEtotSW/totarea
    824             dEtotLW  = dEtotLW/totarea
    825             dEtotsSW = dEtotsSW/totarea
    826             dEtotsLW = dEtotsLW/totarea
     809            dEtotSW  = cpp*SUM(massarea(:,:)*zdtsw(:,:))/totarea
     810            dEtotLW  = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea
     811            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea
     812            dEtotsLW = SUM(fluxsurf_lw(:)*emis(:)*area(:)-zplanck(:))/totarea
     813            dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
     814            dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
    827815            print*,'---------------------------------------------------------------'
    828816            print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
    … …  
    917905         ! test energy conservation
    918906         if(enertest)then
    919             dEtot=0.0
    920             dEtots=0.0
    921             dEzdiff(:,:)=0.
    922             AtmToSurf_TurbFlux=0.
    923             dEdiffs(:)=0.
    924             dEdiff(:)=0.
    925                    
     907            dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
    926908            do ig = 1, ngrid
    927                do l = 1, nlayer
    928                   masse = (pplev(ig,l) - pplev(ig,l+1))/g
    929                   dEzdiff (ig,l)= cpp*masse*(zdtdif(ig,l))
    930                   dEdiff(ig)=dEdiff(ig)+dEzdiff (ig,l)
    931                enddo
     909               dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
    932910               dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
    933                dEdiff(ig)=dEdiff(ig)+ sensibFlux(ig)! subtract flux to the ground
    934                dEtot = dEtot + dEdiff(ig)*area(ig)
    935                dEdiffs(ig)=capcal(ig)*zdtsdif(ig)-zflubid(ig)-sensibFlux(ig)
    936                dEtots = dEtots + dEdiffs(ig)*area(ig)
    937                AtmToSurf_TurbFlux=AtmToSurf_TurbFlux+sensibFlux(ig)*area(ig)
    938911            enddo
    939             dEtot  = dEtot/totarea
    940             dEtots = dEtots/totarea
    941             AtmToSurf_TurbFlux=AtmToSurf_TurbFlux/totarea
     912            dEtot = SUM(dEdiff(:)*area(:))/totarea
     913            dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
     914            dEtots = SUM(dEdiffs(:)*area(:))/totarea
     915            AtmToSurf_TurbFlux=SUM(sensibFlux(:)*area(:))/totarea
    942916            if (UseTurbDiff) then
    943917               print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
    … …  
    954928         !-------------------------
    955929
    956 
    957930         !-------------------------
    958931         ! test water conservation
    959932         if(watertest.and.water)then
    960             dWtot=0.0
    961             dWtots=0.0
    962             nconsMAX=0.0
     933            dWtot = SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap))*ptimestep/totarea
     934            dWtots =  SUM(zdqsdif(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
    963935            do ig = 1, ngrid
    964 
    965                vdifcncons(ig)=0.0
    966                do l = 1, nlayer
    967                   masse = (pplev(ig,l) - pplev(ig,l+1))/g
    968 
    969                   iq    = igcm_h2o_vap
    970                   dWtot = dWtot + masse*zdqdif(ig,l,iq)*ptimestep*area(ig)
    971                   vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
    972 
    973                   iq    = igcm_h2o_ice
    974                   dWtot = dWtot + masse*zdqdif(ig,l,iq)*ptimestep*area(ig)
    975                   vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
    976                
    977                enddo
    978 
    979                iq     = igcm_h2o_vap
    980                dWtots = dWtots + zdqsdif(ig,iq)*ptimestep*area(ig)
    981                vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
    982 
    983                iq     = igcm_h2o_ice
    984                dWtots = dWtots + zdqsdif(ig,iq)*ptimestep*area(ig)
    985                vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)
    986 
    987                if(vdifcncons(ig).gt.nconsMAX)then
    988                   nconsMAX=vdifcncons(ig)
    989                endif
    990 
    991             enddo
    992 
    993             dWtot  = dWtot/totarea
    994             dWtots = dWtots/totarea
     936               vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
     937            Enddo
     938            dWtot = dWtot + SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice))*ptimestep/totarea
     939            dWtots = dWtots + SUM(zdqsdif(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
     940            do ig = 1, ngrid
     941               vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
     942            Enddo           
     943            nconsMAX=MAXVAL(vdifcncons(:))
     944
    995945            print*,'---------------------------------------------------------------'
    996946            print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
    … …  
    998948            print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
    999949            print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
    1000 
    1001950
    1002951         endif
    … …  
    10611010         ! test energy conservation
    10621011         if(enertest)then
    1063             dEtot=0.0
    1064             do ig = 1, ngrid
    1065                do l = 1, nlayer
    1066                   masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1067                   dEtot = dEtot + cpp*masse*zdtadj(ig,l)*area(ig)
    1068                enddo
    1069             enddo
    1070             dEtot=dEtot/totarea
     1012            dEtot=cpp*SUM(massarea(:,:)*zdtadj(:,:))/totarea
    10711013          print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
    10721014         endif
    … …  
    10761018         ! test water conservation
    10771019         if(watertest)then
    1078             dWtot=0.0
     1020            dWtot = SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap))*ptimestep/totarea
    10791021            do ig = 1, ngrid
    1080                cadjncons(ig)=0.0
    1081                do l = 1, nlayer
    1082                   masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1083                  
    1084                   iq    = igcm_h2o_vap
    1085                   dWtot = dWtot + masse*zdqadj(ig,l,iq)*ptimestep*area(ig)
    1086                   cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
    1087                  
    1088                   iq    = igcm_h2o_ice
    1089                   dWtot = dWtot + masse*zdqadj(ig,l,iq)*ptimestep*area(ig)
    1090                   cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
    1091                  
    1092                enddo
    1093 
    1094             enddo
    1095             dWtot=dWtot/totarea
     1022               cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
     1023            Enddo
     1024            dWtot = dWtot + SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice))*ptimestep/totarea
     1025            do ig = 1, ngrid
     1026               cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
     1027            Enddo           
     1028            nconsMAX=MAXVAL(cadjncons(:))
     1029
    10961030            print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
     1031            print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
    10971032         endif
    10981033         !-------------------------
    … …  
    11411076         ! test energy conservation
    11421077         if(enertest)then
    1143             dEtot=0.0
    1144             dEtots=0.0
    1145             do ig = 1, ngrid
    1146                do l = 1, nlayer
    1147                   masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1148                   dEtot = dEtot + cpp*masse*zdtc(ig,l)*area(ig)
    1149                enddo
    1150                dEtots = dEtots + capcal(ig)*zdtsurfc(ig)*area(ig)
    1151             enddo
    1152             dEtot=dEtot/totarea
    1153             dEtots=dEtots/totarea
     1078            dEtot = cpp*SUM(massarea(:,:)*zdtc(:,:))/totarea
     1079            dEtots = SUM(capcal(:)*zdtsurfc(:)*area(:))/totarea
    11541080            print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
    11551081            print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
    … …  
    11991125               ! test energy conservation
    12001126               if(enertest)then
    1201                   dEtot=0.0
    1202                   madjdE(:)=0.0
     1127                  dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
    12031128                  do ig = 1, ngrid
    1204                      do l = 1, nlayer
    1205                         masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1206                         dEtot = dEtot + cpp*masse*(dtmoist(ig,l)+dtevap(ig,l))*area(ig)
    1207                         madjdE(ig) = madjdE(ig) + cpp*masse*(dtmoist(ig,l)+dtevap(ig,l))
    1208                      enddo
     1129                     madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
    12091130                  enddo
    1210                   dEtot=dEtot/totarea
    12111131                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
    1212                endif
    1213                !-------------------------
    1214 
    1215                !-------------------------
    1216                ! test water conservation
    1217                if(watertest)then
    1218                   dWtot=0.0
    1219                   do ig = 1, ngrid
    1220                      do iq = 1 , nq
    1221                         do l = 1, nlayer
    1222                            masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1223                            dWtot = dWtot + masse*dqmoist(ig,l,igcm_h2o_vap)*area(ig)*ptimestep
    1224                            dWtot = dWtot + masse*dqmoist(ig,l,igcm_h2o_ice)*area(ig)*ptimestep
    1225                         enddo
    1226                      enddo
    1227                   enddo
    1228                   dWtot=dWtot/totarea
     1132                 
     1133                ! test energy conservation
     1134                  dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
     1135                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
    12291136                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
    12301137               endif
    … …  
    12641171               ! test energy conservation
    12651172               if(enertest)then
    1266                   dEtot=0.0
    1267                   lscaledE(:)=0.0
     1173                  dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
    12681174                  do ig = 1, ngrid
    1269                      do l = 1, nlayer
    1270                         masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1271                         dEtot = dEtot + cpp*masse*(dtlscale(ig,l)+dtevap(ig,l))*area(ig)
    1272                         lscaledE(ig) = lscaledE(ig) + cpp*masse*(dtlscale(ig,l)+dtevap(ig,l))
    1273                      enddo
     1175                     madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
    12741176                  enddo
    1275                   dEtot=dEtot/totarea
    1276                   print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
    1277                endif
    1278                !-------------------------
    1279 
    1280                !-------------------------
    1281                ! test water conservation
    1282                if(watertest)then
    1283                   dWtot=0.0
    1284                   do ig = 1, ngrid
    1285                      do iq = 1 , nq
    1286                         do l = 1, nlayer
    1287                            masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1288                            dWtot = dWtot + masse*dqvaplscale(ig,l)*area(ig)*ptimestep
    1289                            dWtot = dWtot + masse*dqcldlscale(ig,l)*area(ig)*ptimestep
    1290                         enddo
    1291                      enddo
    1292                   enddo
    1293                   dWtot=dWtot/totarea
    1294                   print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
     1177                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
     1178                 
     1179                ! test energy conservation
     1180                  dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
     1181                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
     1182                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
    12951183               endif
    12961184               !-------------------------
    … …  
    13361224
    13371225
    1338               !-------------------------
    1339               ! test energy conservation
    1340               if(enertest)then
    1341                  dEtot=0.0
    1342                  do ig = 1, ngrid
    1343                     do l = 1, nlayer
    1344                        masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1345                        dEtot = dEtot + cpp*masse*zdtrain(ig,l)*area(ig)
    1346                     enddo
    1347                  enddo
    1348                  dEtot=dEtot/totarea
     1226
     1227               !-------------------------
     1228               ! test energy conservation
     1229               if(enertest)then
     1230                  dEtot=cpp*SUM(massarea(:,:)*zdtrain(:,:))/totarea
    13491231                 print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
    1350 
    1351                  dEtot=0.0
    1352                  do ig = 1, ngrid
    1353                     do l = 1, nlayer
    1354                        masse = (pplev(ig,l) - pplev(ig,l+1))/g
    1355                        dItot = dItot + masse*zdqrain(ig,l,igcm_h2o_ice)*area(ig)
    1356                        dVtot = dVtot + masse*zdqrain(ig,l,igcm_h2o_vap)*area(ig)
    1357                     enddo
    1358                     dItot = dItot + zdqssnow(ig)*area(ig)
    1359                     dVtot = dVtot + zdqsrain(ig)*area(ig)
    1360                  enddo
    1361                  dEtot=(dItot*RLVTT/cpp + dVtot*RLVTT/cpp)/totarea
    1362                  print*,'In rain dItot =',dItot*RLVTT/(cpp*totarea),' W m-2'
    1363                  print*,'In rain dVtot =',dVtot*RLVTT/(cpp*totarea),' W m-2'
     1232                  dItot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))/totarea*RLVTT/cpp
     1233                  dItot = dItot + SUM(area(:)*zdqssnow(:))/totarea*RLVTT/cpp
     1234                  dVtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
     1235                  dVtot = dItot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
     1236                  dEtot = dItot + dVtot
     1237                 print*,'In rain dItot =',dItot,' W m-2'
     1238                 print*,'In rain dVtot =',dVtot,' W m-2'
    13641239                 print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
    1365               endif
    1366               !-------------------------
    1367 
    1368               !-------------------------
    1369               ! test water conservation
    1370               if(watertest)then
    1371                  dWtot=0.0
    1372                  dWtots=0.0
    1373                  do ig = 1, ngrid
    1374                     !do iq = 1 , nq
    1375                        do l = 1, nlayer
    1376                           masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
    1377                           dWtot = dWtot + masse*zdqrain(ig,l,igcm_h2o_vap)*area(ig)*ptimestep
    1378                           dWtot = dWtot + masse*zdqrain(ig,l,igcm_h2o_ice)*area(ig)*ptimestep
    1379                        enddo
    1380                     !enddo
    1381                     dWtots = dWtots + (zdqsrain(ig)+zdqssnow(ig))*area(ig)*ptimestep
    1382                  enddo
    1383                  dWtot=dWtot/totarea
    1384                  dWtots=dWtots/totarea
     1240
     1241               ! test water conservation
     1242                  dWtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
     1243                  dWtot = dWtot + SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))*ptimestep/totarea
     1244                  dWtots =  SUM((zdqsrain(:)+zdqssnow(:))*area(:))*ptimestep/totarea
    13851245                 print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
    13861246                 print*,'In rain surface water change            =',dWtots,' kg m-2'
    … …  
    14061266           ! find qtot
    14071267           if(watertest)then
    1408               dWtot=0.0
    1409               dWtots=0.0
    14101268              iq=3
    1411               do ig = 1, ngrid
    1412                  do l = 1, nlayer
    1413                     masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
    1414                     dWtot  = dWtot  + masse*pq(ig,l,iq)*area(ig)*ptimestep
    1415                     dWtots = dWtots + masse*pdq(ig,l,iq)*area(ig)*ptimestep
    1416                  enddo
    1417               enddo
    1418               dWtot=dWtot/totarea
    1419               dWtots=dWtots/totarea
     1269              dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
     1270              dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
    14201271              print*,'Before sedim pq  =',dWtot,' kg m-2'
    14211272              print*,'Before sedim pdq =',dWtots,' kg m-2'
    … …  
    14291280           ! find qtot
    14301281           if(watertest)then
    1431               dWtot=0.0
    1432               dWtots=0.0
    14331282              iq=3
    1434               do ig = 1, ngrid
    1435                  do l = 1, nlayer
    1436                     masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
    1437                     dWtot  = dWtot  + masse*pq(ig,l,iq)*area(ig)*ptimestep
    1438                     dWtots = dWtots + masse*pdq(ig,l,iq)*area(ig)*ptimestep
    1439                  enddo
    1440               enddo
    1441               dWtot=dWtot/totarea
    1442               dWtots=dWtots/totarea
     1283              dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
     1284              dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
    14431285              print*,'After sedim pq  =',dWtot,' kg m-2'
    14441286              print*,'After sedim pdq =',dWtots,' kg m-2'
    … …  
    14611303           ! test water conservation
    14621304           if(watertest)then
    1463               dWtot=0.0
    1464               dWtots=0.0
    1465               do iq=1,nq
    1466               do ig = 1, ngrid
    1467                  do l = 1, nlayer
    1468                     masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
    1469                     dWtot = dWtot + masse*zdqsed(ig,l,iq)*area(ig)*ptimestep
    1470                  enddo
    1471                  dWtots   = dWtots + zdqssed(ig,iq)*area(ig)*ptimestep
    1472               enddo
    1473               enddo
    1474               dWtot=dWtot/totarea
    1475               dWtots=dWtots/totarea
     1305              dWtot = SUM(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice)))*ptimestep/totarea
     1306              dWtots = SUM((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:))*ptimestep/totarea
    14761307              print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
    14771308              print*,'In sedim surface ice change             =',dWtots,' kg m-2'
    … …  
    15071338            ! test energy conservation
    15081339            if(enertest)then
    1509                dEtot=0.0
    1510                do ig = 1, ngrid
    1511                   dEtots = dEtots + capcal(ig)*zdtsurf_hyd(ig)*area(ig)
    1512                enddo
    1513                dEtot=dEtot/totarea
    1514                print*,'In hydrol atmospheric energy change     =',dEtot,' W m-2'
     1340               dEtots = SUM(area(:)*capcal(:)*zdtsurf_hyd(:))/totarea
     1341               print*,'In hydrol surface energy change     =',dEtots,' W m-2'
    15151342            endif
    15161343            !-------------------------
    … …  
    15191346            ! test water conservation
    15201347            if(watertest)then
    1521                dWtots=0.0
    1522                do ig = 1, ngrid
    1523                   dWtots = dWtots + dqs_hyd(ig,igcm_h2o_ice)*area(ig)*ptimestep
    1524                enddo
    1525                dWtots=dWtots/totarea
     1348               dWtots =  SUM(dqs_hyd(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
    15261349               print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
    1527                dWtots=0.0
    1528                do ig = 1, ngrid
    1529                   dWtots = dWtots + dqs_hyd(ig,igcm_h2o_vap)*area(ig)*ptimestep
    1530                enddo
    1531                dWtots=dWtots/totarea
     1350               dWtots =  SUM(dqs_hyd(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
    15321351               print*,'In hydrol surface water change          =',dWtots,' kg m-2'
    15331352               print*,'---------------------------------------------------------------'
    … …  
    16481467!     ---------------------------------------------------------
    16491468
    1650       Ts1 = 0.0
    1651       Ts2 = 99999.9
    1652       Ts3 = 0.0
    1653       TsS = 0.0 ! mean temperature at bottom soil layer
    1654       do ig=1,ngrid
    1655          Ts1 = Ts1 + area(ig)*tsurf(ig)
    1656          Ts2 = min(Ts2,tsurf(ig)) 
    1657          Ts3 = max(Ts3,tsurf(ig))
    1658          TsS = TsS + area(ig)*tsoil(ig,nsoilmx)
    1659       end do
    1660       Ts1=Ts1/totarea
    1661       TsS=TsS/totarea
     1469      Ts1 = SUM(area(:)*tsurf(:))/totarea
     1470      Ts2 = MINVAL(tsurf(:))
     1471      Ts3 = MAXVAL(tsurf(:))
    16621472      if(callsoil)then
     1473         TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
    16631474         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
    16641475         print*,Ts1,Ts2,Ts3,TsS
    … …  
    16741485      if(corrk)then
    16751486
    1676          ISR = 0.0
    1677          ASR = 0.0
    1678          OLR = 0.0
    1679          GND = 0.0
    1680          DYN = 0.0
    1681          do ig=1,ngrid
    1682             ISR = ISR + area(ig)*fluxtop_dn(ig)
    1683             ASR = ASR + area(ig)*fluxabs_sw(ig)
    1684             OLR = OLR + area(ig)*fluxtop_lw(ig)
    1685             GND = GND + area(ig)*fluxgrd(ig)
    1686             DYN = DYN + area(ig)*fluxdyn(ig)
    1687            
     1487         ISR = SUM(area(:)*fluxtop_dn(:))/totarea
     1488         ASR = SUM(area(:)*fluxabs_sw(:))/totarea
     1489         OLR = SUM(area(:)*fluxtop_lw(:))/totarea
     1490         GND = SUM(area(:)*fluxgrd(:))/totarea
     1491         DYN = SUM(area(:)*fluxdyn(:))/totarea
     1492         do ig=1,ngrid             
    16881493            if(fluxtop_dn(ig).lt.0.0)then
    16891494               print*,'fluxtop_dn has gone crazy'
    … …  
    17021507
    17031508         print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
    1704          print*, ISR/totarea,ASR/totarea,OLR/totarea,GND/totarea,DYN/totarea
     1509         print*, ISR,ASR,OLR,GND,DYN
    17051510         
    17061511         if(enertest)then
    1707             print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
    1708             print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR/totarea,' W m-2'
    1709             print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR/totarea,' W m-2'
     1512            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
     1513            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
     1514            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
    17101515         endif
    17111516
    … …  
    17141519               ! to record global radiative balance
    17151520               open(92,file="rad_bal.out",form='formatted',position='append')
    1716                write(92,*) zday,ISR/totarea,ASR/totarea,OLR/totarea
     1521               write(92,*) zday,ISR,ASR,OLR
    17171522               close(92)
    17181523               open(93,file="tem_bal.out",form='formatted',position='append')
    … …  
    17231528
    17241529      endif
     1530
    17251531
    17261532!     ------------------------------------------------------------------
    … …  
    17771583      if(water)then
    17781584
    1779          icesrf = 0.0
    1780          liqsrf = 0.0
    1781          icecol = 0.0
    1782          vapcol = 0.0
    1783 
    1784          h2otot = 0.0
    1785          do ig=1,ngrid
    1786 
    1787             icesrf = icesrf + area(ig)*qsurf_hist(ig,igcm_h2o_ice)
    1788             liqsrf = liqsrf + area(ig)*qsurf_hist(ig,igcm_h2o_vap)
    1789             icecol = icecol + area(ig)*qcol(ig,igcm_h2o_ice)
    1790             vapcol = vapcol + area(ig)*qcol(ig,igcm_h2o_vap)
    1791 
    1792             h2otot = h2otot + area(ig)*                       &
    1793             (qcol(ig,igcm_h2o_ice)+qcol(ig,igcm_h2o_vap)      &
    1794             +qsurf_hist(ig,igcm_h2o_ice)+qsurf_hist(ig,igcm_h2o_vap))
    1795          end do
    1796 
    1797          print*,' Total water amount [kg m^-2]: ',h2otot/totarea
     1585         icesrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_ice))/totarea
     1586         liqsrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_vap))/totarea
     1587         icecol = SUM(area(:)*qcol(:,igcm_h2o_ice))/totarea
     1588         vapcol = SUM(area(:)*qcol(:,igcm_h2o_vap))/totarea
     1589
     1590         h2otot = icesrf + liqsrf + icecol + vapcol
     1591
     1592         print*,' Total water amount [kg m^-2]: ',h2otot
    17981593         print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
    1799          print*, icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
     1594         print*, icesrf,liqsrf,icecol,vapcol
    18001595
    18011596         if(meanOLR)then
    … …  
    18031598               ! to record global water balance
    18041599               open(98,file="h2o_bal.out",form='formatted',position='append')
    1805                write(98,*) zday,icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
     1600               write(98,*) zday,icesrf,liqsrf,icecol,vapcol
    18061601               close(98)
    18071602            endif
Note: See TracChangeset for help on using the changeset viewer.