Merge pull request #4 from jinyun1tang/jrgrant

Jrgrant new code for maize

examples/sample_run/s022852

@@ -10,7 +10,7 @@
 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00
 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00
 5.77,5.77,5.77,5.77,5.80,5.80,5.80,5.80,4.35,4.35,4.35,4.35,4.35,4.35,4.35,0.00,0.00,0.00,0.00
-0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,5.00,5.00,5.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00
+10.00,10.00,10.00,10.00,10.00,10.00,10.00,10.00,5.00,5.00,5.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00
 3.00,3.00,3.00,3.00,3.00,3.00,3.00,3.00,3.00,3.00,3.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00
 3.87,2.87,2.6,2.6,2.00,2.00,2.00,2.00,3.00,0.50,0.50,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00
 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00

f77src/Makefile

@@ -56,7 +56,7 @@ CMD =		ecosys.x
 # WARNING:  SIMULTANEOUSLY PROFILING AND FLOWTRACING IS NOT RECOMMENDED
 #FFLAGS = -axP -xW -ipo -O3 -r8 -i4 -align dcommons -cpp -save
 #intel setup
-FFLAGS = -O2 -mp1 -r8 -i4 -align dcommons -cpp -auto-scalar
+FFLAGS = -O2 -mp1 -r8 -i4 -align dcommons -cpp -auto-scalar -fimf-arch-consistency=true
 #gfortran
 #FFLAGS = -O2
 

f77src/blk13a.h

@@ -19,7 +19,7 @@
      1,ZNO2SH(JZ,JY,JX),Z2GSH(JZ,JY,JX),Z2OSH(JZ,JY,JX),TRC0(JY,JX)
      2,ZNO2BH(JZ,JY,JX),RC0(0:5,JY,JX),RA0(0:5,JY,JX),TRA0(JY,JX)
      3,TOQCK(0:JZ,JY,JX),ORGN(0:JZ,JY,JX),ORGR(0:JZ,JY,JX)
-     4,ZNFNI(0:JZ,JY,JX),ZNFNG(0:JZ,JY,JX),ZNFN0(0:JZ,JY,JX)
+     4,ZNFNI(0:JZ,JY,JX),ZNFN0(0:JZ,JY,JX)
      5,ZNHUI(0:JZ,JY,JX),ZNHU0(0:JZ,JY,JX),H1PO4(0:JZ,JY,JX)
      6,H1POB(0:JZ,JY,JX),H1PO4H(JZ,JY,JX),H1POBH(JZ,JY,JX)
 

f77src/blk8b.h

@@ -6,6 +6,6 @@
      6,YDPTH(JZ,JY,JX),POROQ(0:JZ,JY,JX),TFND(0:JZ,JY,JX),VOLXA(JY,JX)
      7,PSIMS(JY,JX),PSIMX(JY,JX),PSIMN(JY,JX),PSISD(JY,JX),PSIMD(JY,JX)
      8,SAND(JZ,JY,JX),SILT(JZ,JY,JX),CLAY(JZ,JY,JX),CDPTHZ(0:JZ,JY,JX)
-     9,DPTHZ(JZ,JY,JX),AREA(3,0:JZ,JY,JX),DISP(3,JD,JV,JH),OXKM
+     9,DPTHZ(JZ,JY,JX),AREA(3,0:JZ,JY,JX),DISP(3,JD,JV,JH),OXKM,PSIHY
      1,OMCI(3,0:4),OMCF(7),OMCA(7)
      2,IUTYP(JY,JX),IXTYP(2,JY,JX),IYTYP(0:2,366,JY,JX)

f77src/grosub.f

@@ -40,7 +40,7 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
       DIMENSION RTNT(2),RLNT(2,JZ),RTSK1(2,JZ,10),RTSK2(2,JZ,10)
      2,RTDPL(10,JZ),FWTR(JZ),FWTB(JP),FRTDP(0:3),RCCX(0:2),RCCQ(0:2)
      3,RCCZ(0:3),RCCY(0:3),EFIRE(2,5:5),WGLFBL(JZ,10,JP,JY,JX)
-     4,WTSHTA(JZ,JY,JX),FLG4Y(0:3),ATRPX(0:1),GVMX(0:1)
+     4,WTSHTA(JZ,JY,JX),FLG4Y(0:3),ATRPX(0:1),GVMX(0:1),RTSK(0:2)
       DIMENSION CH2O3(25),CH2O4(25),CPOOLK(10,JP,JY,JX),FHVSTK(0:25)
      2,FHVSHK(0:25),WFNGR(2,JZ),PSILY(0:2)
       DIMENSION FWOOD(0:1),FWOODN(0:1),FWOODP(0:1)
@@ -71,6 +71,7 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
       DATA RCCY/0.333,0.333,0.333,0.333/
       DATA RCCX/0.250,0.833,0.833/
       DATA RCCQ/0.833,0.833,0.833/
+      DATA RTSK/1.0,1.0,4.0/
       DATA FXRN/0.50,0.05,0.50,0.05/
       DATA FXFB/1.0E-02,1.0E-02,1.0E-05,1.0E-05/
       DATA FPART1/1.00/,FPART2/0.40/
@@ -1373,14 +1374,17 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
 C
       IF(IDAY(1,NB,NZ,NY,NX).NE.0
      2.AND.CCPOLB(NB,NZ,NY,NX).GT.ZERO)THEN
-      CCC=AMIN1(CZPOLB(NB,NZ,NY,NX)/(CZPOLB(NB,NZ,NY,NX)
+      CCC=AMAX1(0.0,AMIN1(1.0
+     1,CZPOLB(NB,NZ,NY,NX)/(CZPOLB(NB,NZ,NY,NX)
      2+CCPOLB(NB,NZ,NY,NX)*CNKI)
      3,CPPOLB(NB,NZ,NY,NX)/(CPPOLB(NB,NZ,NY,NX)
-     4+CCPOLB(NB,NZ,NY,NX)*CPKI))
-      CNC=CCPOLB(NB,NZ,NY,NX)/(CCPOLB(NB,NZ,NY,NX)
-     2+CZPOLB(NB,NZ,NY,NX)/CNKI)
-      CPC=CCPOLB(NB,NZ,NY,NX)/(CCPOLB(NB,NZ,NY,NX)
-     2+CPPOLB(NB,NZ,NY,NX)/CPKI)
+     4+CCPOLB(NB,NZ,NY,NX)*CPKI)))
+      CNC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLB(NB,NZ,NY,NX)/(CCPOLB(NB,NZ,NY,NX)
+     2+CZPOLB(NB,NZ,NY,NX)/CNKI)))
+      CPC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLB(NB,NZ,NY,NX)/(CCPOLB(NB,NZ,NY,NX)
+     2+CPPOLB(NB,NZ,NY,NX)/CPKI)))
       ELSE
       CCC=0.0
       CNC=0.0
@@ -2530,16 +2534,17 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
 C     OF STALK RESERVES
 C
       IF(WTRSVB(NB,NZ,NY,NX).GT.ZEROP(NZ,NY,NX))THEN
-      ZNPG=AMIN1(WTRSBN(NB,NZ,NY,NX)/(WTRSBN(NB,NZ,NY,NX)
+      ZNPGN=WTRSBN(NB,NZ,NY,NX)/(WTRSBN(NB,NZ,NY,NX)
      2+SETN*WTRSVB(NB,NZ,NY,NX))
-     3,WTRSBP(NB,NZ,NY,NX)/(WTRSBP(NB,NZ,NY,NX)
-     3+SETP*WTRSVB(NB,NZ,NY,NX)))
-      ZPGRX=ZPGRM+ZPGRD*AMAX1(0.0,AMIN1(1.0,ZNPG))
+      ZNPGP=WTRSBP(NB,NZ,NY,NX)/(WTRSBP(NB,NZ,NY,NX)
+     3+SETP*WTRSVB(NB,NZ,NY,NX))
+      ZPGRN=ZPGRM+ZPGRD*AMAX1(0.0,AMIN1(1.0,ZNPGN))
+      ZPGRP=ZPGRM+ZPGRD*AMAX1(0.0,AMIN1(1.0,ZNPGP))
       XLOCN=AMIN1(XLOCM*CNGR(NZ,NY,NX)
-     2,AMAX1(0.0,WTRSBN(NB,NZ,NY,NX)*ZPGRX)
+     2,AMAX1(0.0,WTRSBN(NB,NZ,NY,NX)*ZPGRN)
      3,(WTGRB(NB,NZ,NY,NX)+XLOCC)*CNGR(NZ,NY,NX)-WTGRBN(NB,NZ,NY,NX))
       XLOCP=AMIN1(XLOCM*CPGR(NZ,NY,NX)
-     2,AMAX1(0.0,WTRSBP(NB,NZ,NY,NX)*ZPGRX)
+     2,AMAX1(0.0,WTRSBP(NB,NZ,NY,NX)*ZPGRP)
      3,(WTGRB(NB,NZ,NY,NX)+XLOCC)*CPGR(NZ,NY,NX)-WTGRBP(NB,NZ,NY,NX))
       ELSE
       XLOCN=0.0
@@ -3233,11 +3238,15 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
       CPPOLN=1.0
       ENDIF
       IF(CCPOLN.GT.ZERO)THEN 
-      CCC=AMIN1(CZPOLN/(CZPOLN+CCPOLN*CNKI)
-     2,CPPOLN/(CPPOLN+CCPOLN*CPKI))
-      CNC=CCPOLN/(CCPOLN+CZPOLN/CNKI)
-      CPC=CCPOLN/(CCPOLN+CPPOLN/CPKI)
-      CNF=CCPOLN/(CCPOLN+CZPOLN/CNKF)
+      CCC=AMAX1(0.0,AMIN1(1.0
+     1,CZPOLN/(CZPOLN+CCPOLN*CNKI)
+     2,CPPOLN/(CPPOLN+CCPOLN*CPKI)))
+      CNC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLN/(CCPOLN+CZPOLN/CNKI)))
+      CPC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLN/(CCPOLN+CPPOLN/CPKI)))
+      CNF=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLN/(CCPOLN+CZPOLN/CNKF)))
       ELSE
       CCC=0.0
       CNC=0.0
@@ -3369,8 +3378,7 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
       IF(CPOOL(NB,NZ,NY,NX).GT.ZEROP(NZ,NY,NX)
      2.AND.WTLSB(NB,NZ,NY,NX).GT.ZEROL(NZ,NY,NX))THEN
       WTLSB1=WTLSB(NB,NZ,NY,NX)
-      WTNDB1=AMIN1(WTLSB(NB,NZ,NY,NX),AMAX1(FSNKM
-     2*WTLSB(NB,NZ,NY,NX),WTNDB(NB,NZ,NY,NX)))
+      WTNDB1=AMIN1(WTLSB(NB,NZ,NY,NX),WTNDB(NB,NZ,NY,NX))
       WTLSBT=WTLSB1+WTNDB1
       IF(WTLSBT.GT.ZEROP(NZ,NY,NX))THEN
       CPOOLD=(CPOOL(NB,NZ,NY,NX)*WTNDB1
@@ -3480,7 +3488,8 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
       IF(RTDP1(N,NR,NZ,NY,NX).GT.CDPTHZ(L-1,NY,NX))THEN
       IF(RTDP1(N,NR,NZ,NY,NX).LE.CDPTHZ(L,NY,NX))THEN
       RTDPP=RTDP1(N,NR,NZ,NY,NX)+HTSTZ(NZ,NY,NX)
-      RTSK1(N,L,NR)=XRTN1*RRAD1(N,L,NZ,NY,NX)**2/RTDPP
+      RTSK1(N,L,NR)=RTSK(IGTYP(NZ,NY,NX))*XRTN1
+     2*RRAD1(N,L,NZ,NY,NX)**2/RTDPP
       RTNT(N)=RTNT(N)+RTSK1(N,L,NR)
       RLNT(N,L)=RLNT(N,L)+RTSK1(N,L,NR)
       ENDIF
@@ -3633,14 +3642,17 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
 C
       IF(IDAY(1,NB1(NZ,NY,NX),NZ,NY,NX).NE.0
      2.AND.CCPOLR(N,L,NZ,NY,NX).GT.ZERO)THEN
-      CCC=AMIN1(CZPOLR(N,L,NZ,NY,NX)/(CZPOLR(N,L,NZ,NY,NX)
+      CCC=AMAX1(0.0,AMIN1(1.0
+     1,CZPOLR(N,L,NZ,NY,NX)/(CZPOLR(N,L,NZ,NY,NX)
      2+CCPOLR(N,L,NZ,NY,NX)*CNKI)
      3,CPPOLR(N,L,NZ,NY,NX)/(CPPOLR(N,L,NZ,NY,NX)
-     4+CCPOLR(N,L,NZ,NY,NX)*CPKI))
-      CNC=CCPOLR(N,L,NZ,NY,NX)/(CCPOLR(N,L,NZ,NY,NX)
-     2+CZPOLR(N,L,NZ,NY,NX)/CNKI)
-      CPC=CCPOLR(N,L,NZ,NY,NX)/(CCPOLR(N,L,NZ,NY,NX)
-     2+CPPOLR(N,L,NZ,NY,NX)/CPKI)
+     4+CCPOLR(N,L,NZ,NY,NX)*CPKI)))
+      CNC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLR(N,L,NZ,NY,NX)/(CCPOLR(N,L,NZ,NY,NX)
+     2+CZPOLR(N,L,NZ,NY,NX)/CNKI)))
+      CPC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLR(N,L,NZ,NY,NX)/(CCPOLR(N,L,NZ,NY,NX)
+     2+CPPOLR(N,L,NZ,NY,NX)/CPKI)))
       ELSE
       CCC=0.0
       CNC=0.0
@@ -3869,14 +3881,17 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
 C
       IF(IDAY(1,NB1(NZ,NY,NX),NZ,NY,NX).NE.0
      2.AND.CCPOLR(N,L,NZ,NY,NX).GT.ZERO)THEN
-      CCC=AMIN1(CZPOLR(N,L,NZ,NY,NX)/(CZPOLR(N,L,NZ,NY,NX)
+      CCC=AMAX1(0.0,AMIN1(1.0
+     1,CZPOLR(N,L,NZ,NY,NX)/(CZPOLR(N,L,NZ,NY,NX)
      2+CCPOLR(N,L,NZ,NY,NX)*CNKI)
      3,CPPOLR(N,L,NZ,NY,NX)/(CPPOLR(N,L,NZ,NY,NX)
-     4+CCPOLR(N,L,NZ,NY,NX)*CPKI))
-      CNC=CCPOLR(N,L,NZ,NY,NX)/(CCPOLR(N,L,NZ,NY,NX)
-     2+CZPOLR(N,L,NZ,NY,NX)/CNKI)
-      CPC=CCPOLR(N,L,NZ,NY,NX)/(CCPOLR(N,L,NZ,NY,NX)
-     2+CPPOLR(N,L,NZ,NY,NX)/CPKI)
+     4+CCPOLR(N,L,NZ,NY,NX)*CPKI)))
+      CNC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLR(N,L,NZ,NY,NX)/(CCPOLR(N,L,NZ,NY,NX)
+     2+CZPOLR(N,L,NZ,NY,NX)/CNKI)))
+      CPC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLR(N,L,NZ,NY,NX)/(CCPOLR(N,L,NZ,NY,NX)
+     2+CPPOLR(N,L,NZ,NY,NX)/CPKI)))
       ELSE
       CCC=0.0
       CNC=0.0
@@ -4475,11 +4490,15 @@ SUBROUTINE grosub(I,J,NHW,NHE,NVN,NVS)
       CPPOLN=1.0
       ENDIF
       IF(CCPOLN.GT.ZERO)THEN 
-      CCC=AMIN1(CZPOLN/(CZPOLN+CCPOLN*CNKI)
-     2,CPPOLN/(CPPOLN+CCPOLN*CPKI))
-      CNC=CCPOLN/(CCPOLN+CZPOLN/CNKI)
-      CPC=CCPOLN/(CCPOLN+CPPOLN/CPKI)
-      CNF=CCPOLN/(CCPOLN+CZPOLN/CNKF)
+      CCC=AMAX1(0.0,AMIN1(1.0
+     1,CZPOLN/(CZPOLN+CCPOLN*CNKI)
+     2,CPPOLN/(CPPOLN+CCPOLN*CPKI)))
+      CNC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLN/(CCPOLN+CZPOLN/CNKI)))
+      CPC=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLN/(CCPOLN+CPPOLN/CPKI)))
+      CNF=AMAX1(0.0,AMIN1(1.0
+     1,CCPOLN/(CCPOLN+CZPOLN/CNKF)))
       ELSE
       CCC=0.0
       CNC=0.0

f77src/hfunc.f

@@ -279,7 +279,7 @@ SUBROUTINE hfunc(I,J,NHW,NHE,NVN,NVS)
 C
       IF(ISTYP(NZ,NY,NX).EQ.0.AND.IDAY(6,NB,NZ,NY,NX).EQ.0)THEN
       WFNS=AMIN1(1.0,AMAX1(0.0,PSILG(NZ,NY,NX)-PSILM))
-      WFNSP=WFNS**0.25
+      WFNSP=WFNS**0.167
       RNI=RNI*WFNSP
       RLA=RLA*WFNSP
       ENDIF