C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagstats_calc.F,v 1.6 2014/08/25 21:50:02 jmc Exp $ C $Name: $ #include "DIAG_OPTIONS.h" CBOP C !ROUTINE: DIAGSTATS_CALC C !INTERFACE: SUBROUTINE DIAGSTATS_CALC( O statArr, I inpArr, frcArr, scaleFact, power, useFract, I useReg, regMskVal, I nStats,sizI1,sizI2,sizJ1,sizJ2, iRun,jRun, I regMask, arrMask, arrhFac, arrArea, I arrDr, specialVal, exclSpVal, useWeight, I myThid ) C !DESCRIPTION: C Compute statistics for this tile, level, region C !USES: IMPLICIT NONE #include "EEPARAMS.h" #include "SIZE.h" C !INPUT/OUTPUT PARAMETERS: C == Routine Arguments == C statArr :: output statistics array C inpArr :: input field array to process (compute stats & add to statFld) C frcArr :: fraction used for weighted-average diagnostics C scaleFact :: scaling factor C power :: option to fill-in with the field square (power=2) C useFract :: if True, use fraction-weight C useReg :: how to use region-mask: =0 : not used ; C =1 : grid-center location ; =2 : U location ; =3 : V location C regMskVal :: region-mask identificator value C (point i,j belong to region <=> regMask(i,j) = regMskVal) C nStats :: size of output array: statArr C sizI1,sizI2 :: size of inpArr array: 1rst index range (min,max) C sizJ1,sizJ2 :: size of inpArr array: 2nd index range (min,max) C iRun,jRun :: range of 1rst & 2nd index to process C regMask :: regional mask C arrMask :: mask for this input array C arrhFac :: weight factor (horizontally varying) C arrArea :: Area weighting factor C arrDr :: uniform weighting factor C specialVal :: special value in input array (to exclude if exclSpVal=T) C exclSpVal :: if T, exclude "specialVal" in input array C useWeight :: use weight factor "arrhFac" Cc k,bi,bj :: level and tile indices used for weighting (mask,area ...) Cc parsFld :: parser field with characteristics of the diagnostics C myThid :: my Thread Id number INTEGER nStats,sizI1,sizI2,sizJ1,sizJ2 INTEGER iRun, jRun _RL statArr(0:nStats) _RL inpArr (sizI1:sizI2,sizJ1:sizJ2) _RL frcArr (sizI1:sizI2,sizJ1:sizJ2) _RL scaleFact INTEGER power LOGICAL useFract INTEGER useReg _RS regMskVal _RS regMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS arrMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS arrhFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS arrArea(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL arrDr _RL specialVal LOGICAL exclSpVal LOGICAL useWeight INTEGER myThid CEOP C !LOCAL VARIABLES: C i,j :: loop indices INTEGER i, j, n INTEGER im, ix #ifndef TARGET_NEC_SX LOGICAL inside(sNx+1,sNy+1) _RL tmpFld(sNx+1,sNy+1) _RL tmpVol(sNx+1,sNy+1) #else C Extra variables and fields to support vectorization. C This code also uses the intrinsic F90 routines SUM, MINVAL, MAXVAL C and thus will not compile with a F77 compiler. _RL arrMaskL(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL tmpVol (1-OLx:sNx+OLx,1-OLy:sNy+OLy) #endif im = nStats - 1 ix = nStats DO n=0,nStats statArr(n) = 0. ENDDO #ifndef TARGET_NEC_SX C- Apply Scaling Factor and power option to Input Field (-> tmpFld): IF ( power.EQ.2 ) THEN DO j = 1,jRun DO i = 1,iRun tmpFld(i,j) = scaleFact*inpArr(i,j) tmpFld(i,j) = tmpFld(i,j)*tmpFld(i,j) ENDDO ENDDO ELSE DO j = 1,jRun DO i = 1,iRun tmpFld(i,j) = scaleFact*inpArr(i,j) ENDDO ENDDO ENDIF C- Set weight factor "tmpVol" (area and hFac and/or fraction field) C and part of domain (=inside) where to compute stats IF ( useFract .AND. useWeight ) THEN DO j = 1,jRun DO i = 1,iRun inside(i,j) = arrMask(i,j).NE.0. & .AND. arrhFac(i,j).NE.0. & .AND. frcArr(i,j) .NE.0. tmpVol(i,j) = arrArea(i,j)*arrhFac(i,j)*frcArr(i,j) ENDDO ENDDO ELSEIF ( useFract ) THEN DO j = 1,jRun DO i = 1,iRun inside(i,j) = arrMask(i,j).NE.0. & .AND. arrhFac(i,j).NE.0. & .AND. frcArr(i,j) .NE.0. tmpVol(i,j) = arrArea(i,j)*frcArr(i,j) ENDDO ENDDO ELSEIF ( useWeight ) THEN DO j = 1,jRun DO i = 1,iRun inside(i,j) = arrMask(i,j).NE.0. & .AND. arrhFac(i,j).NE.0. tmpVol(i,j) = arrArea(i,j)*arrhFac(i,j) ENDDO ENDDO ELSE DO j = 1,jRun DO i = 1,iRun inside(i,j) = arrMask(i,j).NE.0. & .AND. arrhFac(i,j).NE.0. tmpVol(i,j) = arrArea(i,j) ENDDO ENDDO ENDIF C- Exclude (setting inside=F) Special Value: IF ( exclSpVal ) THEN DO j = 1,jRun DO i = 1,iRun inside(i,j) = inside(i,j) .AND. inpArr(i,j).NE.specialVal ENDDO ENDDO ENDIF C- Account for Region-mask (refine "inside"): IF ( useReg.EQ.1 ) THEN DO j = 1,jRun DO i = 1,iRun inside(i,j) = inside(i,j) .AND. regMask(i,j).EQ.regMskVal ENDDO ENDDO ELSEIF ( useReg.EQ.2 ) THEN DO j = 1,jRun DO i = 1,iRun inside(i,j) = inside(i,j) .AND.( regMask(i,j).EQ.regMskVal & .OR. regMask(i-1,j).EQ.regMskVal ) ENDDO ENDDO ELSEIF ( useReg.EQ.3 ) THEN DO j = 1,jRun DO i = 1,iRun inside(i,j) = inside(i,j) .AND.( regMask(i,j).EQ.regMskVal & .OR. regMask(i,j-1).EQ.regMskVal ) ENDDO ENDDO ENDIF C- Calculate Stats DO j = 1,jRun DO i = 1,iRun IF ( inside(i,j) ) THEN statArr(im) = tmpFld(i,j) statArr(0) = statArr(0) + tmpVol(i,j) statArr(1) = statArr(1) + tmpVol(i,j)*tmpFld(i,j) statArr(2) = statArr(2) + tmpVol(i,j)*tmpFld(i,j)*tmpFld(i,j) ENDIF ENDDO ENDDO statArr(ix) = statArr(im) DO j = 1,jRun DO i = 1,iRun IF ( inside(i,j) ) THEN statArr(im) = MIN(tmpFld(i,j),statArr(im)) statArr(ix) = MAX(tmpFld(i,j),statArr(ix)) ENDIF ENDDO ENDDO statArr(0) = statArr(0)*arrDr statArr(1) = statArr(1)*arrDr statArr(2) = statArr(2)*arrDr #else /* TARGET_NEC_SX defined */ arrMaskL = 0. _d 0 IF ( useFract .AND. exclSpVal ) THEN DO j = 1,jRun DO i = 1,iRun IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0. & .AND. arrhFac(i,j).NE.0. & .AND. inpArr(i,j).NE.specialVal ) & arrMaskL(i,j) = 1. _d 0 ENDDO ENDDO IF ( useWeight ) THEN tmpVol = arrhFac*arrArea*frcArr ELSE tmpVol = arrArea*frcArr ENDIF ELSEIF ( useFract ) THEN DO j = 1,jRun DO i = 1,iRun IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0. & .AND. arrhFac(i,j).NE.0. ) & arrMaskL(i,j) = 1. _d 0 ENDDO ENDDO IF ( useWeight ) THEN tmpVol = arrhFac*arrArea*frcArr ELSE tmpVol = arrArea*frcArr ENDIF ELSEIF ( exclSpVal ) THEN DO j = 1,jRun DO i = 1,iRun IF ( arrMask(i,j).NE.0. & .AND. arrhFac(i,j).NE.0. & .AND. inpArr(i,j).NE.specialVal ) & arrMaskL(i,j) = 1. _d 0 ENDDO ENDDO IF ( useWeight ) THEN tmpVol = arrhFac*arrArea ELSE tmpVol = arrArea ENDIF ELSE DO j = 1,jRun DO i = 1,iRun IF ( arrMask(i,j).NE.0. & .AND. arrhFac(i,j).NE.0. ) & arrMaskL(i,j) = 1. _d 0 ENDDO ENDDO IF ( useWeight ) THEN tmpVol = arrhFac*arrArea ELSE tmpVol = arrArea ENDIF ENDIF C- Account for Region-mask: IF ( useReg.EQ.1 ) THEN DO j = 1,jRun DO i = 1,iRun IF ( regMask(i,j).NE.regMskVal ) arrMaskL(i,j) = 0. _d 0 ENDDO ENDDO ELSEIF ( useReg.EQ.2 ) THEN DO j = 1,jRun DO i = 1,iRun IF ( regMask(i,j).NE.regMskVal .AND. & regMask(i-1,j).NE.regMskVal ) arrMaskL(i,j) = 0. _d 0 ENDDO ENDDO ELSEIF ( useReg.EQ.3 ) THEN DO j = 1,jRun DO i = 1,iRun IF ( regMask(i,j).NE.regMskVal .AND. & regMask(i,j-1).NE.regMskVal ) arrMaskL(i,j) = 0. _d 0 ENDDO ENDDO ENDIF C inpArr can be undefined/non-initialised in overlaps, so we need C to clean this fields first by copying the defined range to tmpFld tmpFld = 0. _d 0 DO j = 1,jRun DO i = 1,iRun tmpFld(i,j) = inpArr(i,j)*scaleFact ENDDO ENDDO IF ( power.EQ.2) THEN tmpFld = tmpFld*tmpFld ENDIF C sum up the volume tmpVol = tmpVol*arrMaskL statArr(0) = SUM(tmpVol)*arrDr C compute and sum up volume*field tmpVol = tmpVol*tmpFld statArr(1) = SUM(tmpVol)*arrDr C compute and sum up volume*field**2 tmpVol = tmpVol*tmpFld statArr(2) = SUM(tmpVol)*arrDr statArr(im) = MINVAL(tmpFld, MASK = arrMaskL0.) statArr(ix) = MAXVAL(tmpFld, MASK = arrMaskL0.) #endif /* TARGET_NEC_SX */ RETURN END