forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
mon_calc_stats_rs.F
161 lines (144 loc) · 4.91 KB
/
mon_calc_stats_rs.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_calc_stats_rs.F,v 1.2 2016/07/25 21:48:56 jmc Exp $
C $Name: $
#include "MONITOR_OPTIONS.h"
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: MON_CALC_STATS_RS
C !INTERFACE:
SUBROUTINE MON_CALC_STATS_RS(
I myNr, arr, arrhFac, arrMask, arrArea, arrDr,
O theMin, theMax, theMean, theSD, theDel2, theVol,
I myThid )
C Calculate statistics of global array ``\_RS arr''.
C account for volume and mask
C !USES:
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
C !INPUT/OUTPUT PARAMETERS:
INTEGER myNr
_RS arr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,myNr,nSx,nSy)
_RS arrhFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,myNr,nSx,nSy)
_RS arrMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RS arrArea(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RS arrDr(myNr)
_RL theMin, theMax, theMean, theSD, theDel2, theVol
INTEGER myThid
CEOP
C !LOCAL VARIABLES:
INTEGER bi,bj,i,j,k
LOGICAL noPnts
_RL tmpVal
_RL tmpMask
_RL tmpVol
_RL ddx, ddy
_RL theVar
_RL theNbPt
_RL tileMean(nSx,nSy)
_RL tileVar (nSx,nSy)
_RL tileSD (nSx,nSy)
_RL tileDel2(nSx,nSy)
_RL tileVol (nSx,nSy)
_RL tileNbPt(nSx,nSy)
theMin = 0.
theMax = 0.
theMean= 0.
theSD = 0.
theVar = 0.
theDel2= 0.
theVol = 0.
theNbPt= 0.
noPnts = .TRUE.
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
tileNbPt(bi,bj) = 0.
tileDel2(bi,bj) = 0.
tileVol (bi,bj) = 0.
tileMean(bi,bj) = 0.
tileVar (bi,bj) = 0.
DO k=1,myNr
DO j=1,sNy
DO i=1,sNx
tmpVal = arr(i,j,k,bi,bj)
tmpMask = arrMask(i,j,bi,bj)*arrhFac(i,j,k,bi,bj)
IF ( tmpMask.GT.0. _d 0 .AND. noPnts ) THEN
theMin=tmpVal
theMax=tmpVal
noPnts=.FALSE.
ENDIF
IF ( tmpMask.GT.0. _d 0 ) THEN
theMin = MIN(theMin,tmpVal)
theMax = MAX(theMax,tmpVal)
C-- like old code (but using hFac instead of mask): identical if no partial cell
c tileDel2(bi,bj) = tileDel2(bi,bj)
c & + 0.25*ABS(
c & (arr(i+1,j,k,bi,bj)-tmpVal)*arrhFac(i+1,j,k,bi,bj)
c & +(arr(i-1,j,k,bi,bj)-tmpVal)*arrhFac(i-1,j,k,bi,bj)
c & +(arr(i,j+1,k,bi,bj)-tmpVal)*arrhFac(i,j+1,k,bi,bj)
c & +(arr(i,j-1,k,bi,bj)-tmpVal)*arrhFac(i,j-1,k,bi,bj)
c & )
C-- New form:
ddx = arrhFac(i+1,j,k,bi,bj)*arrhFac(i-1,j,k,bi,bj)
IF ( ddx.GT.0. _d 0 ) THEN
ddx = (arr(i+1,j,k,bi,bj)-tmpVal)
& + (arr(i-1,j,k,bi,bj)-tmpVal)
ENDIF
ddy = arrhFac(i,j+1,k,bi,bj)*arrhFac(i,j-1,k,bi,bj)
IF ( ddy.GT.0. _d 0 ) THEN
ddy = (arr(i,j+1,k,bi,bj)-tmpVal)
& + (arr(i,j-1,k,bi,bj)-tmpVal)
ENDIF
tileDel2(bi,bj) = tileDel2(bi,bj) + ddx*ddx + ddy*ddy
tileNbPt(bi,bj) = tileNbPt(bi,bj) + oneRL
tmpVol = arrArea(i,j,bi,bj)*arrDr(k)*tmpMask
tileVol (bi,bj) = tileVol (bi,bj) + tmpVol
tileMean(bi,bj) = tileMean(bi,bj) + tmpVol*tmpVal
tileVar (bi,bj) = tileVar (bi,bj) + tmpVol*tmpVal*tmpVal
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
CALL GLOBAL_SUM_TILE_RL( tileNbPt, theNbPt, myThid )
CALL GLOBAL_SUM_TILE_RL( tileDel2, theDel2, myThid )
CALL GLOBAL_SUM_TILE_RL( tileVol , theVol , myThid )
CALL GLOBAL_SUM_TILE_RL( tileMean, theMean, myThid )
c CALL GLOBAL_SUM_TILE_RL( tileVar , theVar , myThid )
IF ( theNbPt.GT.zeroRL ) THEN
c theDel2 = theDel2/theNbPt
theDel2 = SQRT(theDel2)/theNbPt
ENDIF
IF ( theVol.GT.0. _d 0 ) THEN
theMean= theMean/theVol
theVar = theVar/theVol
IF ( noPnts ) theMin = theMean
theMin = -theMin
_GLOBAL_MAX_RL(theMin,myThid)
theMin = -theMin
IF ( noPnts ) theMax = theMean
_GLOBAL_MAX_RL(theMax,myThid)
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
tileSD(bi,bj)=0.
DO k=1,myNr
DO j=1,sNy
DO i=1,sNx
tmpVal=arr(i,j,k,bi,bj)
tmpMask = arrMask(i,j,bi,bj)*arrhFac(i,j,k,bi,bj)
IF ( tmpMask.GT.0. _d 0 ) THEN
tmpVol = arrArea(i,j,bi,bj)*arrDr(k)*tmpMask
tileSD(bi,bj) = tileSD(bi,bj)
& + tmpVol*(tmpVal-theMean)*(tmpVal-theMean)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
CALL GLOBAL_SUM_TILE_RL( tileSD, theSD, myThid )
theSD = SQRT(theSD/theVol)
c theSD = SQRT(theVar-theMean*theMean)
ENDIF
RETURN
END