pkg/monitor/mon_surfcor.F

C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_surfcor.F,v 1.4 2003/12/24 00:29:17 jmc Exp $
C $Name:  $

#include "MONITOR_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: MON_SURFCOR

C     !INTERFACE:
      SUBROUTINE MON_SURFCOR(
     I     myThid )

C     !DESCRIPTION:
C     Compute and write area-mean surface expansion term (also called
C     ``surface correction'' with Linear FS).
C
C     Diagnose mean surface expansion term
C     \begin{equation}
C       \mbox{with r coordinate} = (\mbox{w surf})(\mbox{Tracer})
C     \end{equation}
C     \begin{equation}
C       \mbox{units} = (\mbox{W units})(\mbox{Tracer units})
C       \ ; \ \mbox{+ = out} 
C     \end{equation}
C     \begin{equation}
C       \mbox{with r* coord}) = \frac{d\eta}{dt} \frac{dz}{H}
C       (\mbox{Tracer})
C     \end{equation}
C     
C     Atmosphere: convert surf.cor(Theta) to surface heating, 
C     \begin{equation}
C       \mbox{units} = \frac{W}{m^2}, \mbox{+ = out}
C     \end{equation}
C     compute mean conversion term Temp -> PE , units= W/m2, 
C     + = decreasing heat content, increasing PE

C     !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "SURFACE.h"
#include "GRID.h"
#include "MONITOR.h"

C     !INPUT PARAMETERS:
      INTEGER myThid
CEOP

C     !LOCAL VARIABLES:
      INTEGER i,j,k,ks,bi,bj 
      _RL theArea, wT_Mean, wS_Mean, tmp_wS_M, wT_Heat
      _RL vT_Mean, vS_Mean, vT_Heat, theta2PE
      _RL tmpVol, tmpVal, conv_th2Heat, ddPI

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       theArea = 0.
       theta2PE = 0.
       wT_Mean = 0.
       wS_Mean = 0.
       wT_Heat = 0.
       vT_Mean = 0.
       vS_Mean = 0.
       vT_Heat = 0.
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid) 
C-- Compute surface "expansion" term & do the integral
          tmp_wS_M = wS_Mean
          DO j=1,sNy
           DO i=1,sNx
            ks = ksurfC(i,j,bi,bj)
            IF (ks.LE.Nr) THEN
             theArea = theArea + rA(i,j,bi,bj) 
             tmpVal =
     &          rA(i,j,bi,bj)*wVel(i,j,ks,bi,bj)*theta(i,j,ks,bi,bj)
             wT_Mean = wT_Mean + tmpVal
             wS_Mean = wS_Mean
     &        + rA(i,j,bi,bj)*wVel(i,j,ks,bi,bj)* salt(i,j,ks,bi,bj)
C-  Atmos in Pot.Temp => convert Omega*Theta to heat flux :
             IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
               wT_Heat = wT_Heat 
     &                 + tmpVal*atm_cp*((rC(ks)/atm_po)**atm_kappa) 
             ENDIF
            ENDIF
           ENDDO
          ENDDO
#ifdef ALLOW_AIM
          IF ( useAIM ) THEN
           DO j=1,sNy
            DO i=1,sNx
             ks = ksurfC(i,j,bi,bj)
             IF (ks.LE.Nr) THEN 
              tmpVal = salt(i,j,ks,bi,bj)
     &               + salt(i,j,Nr,bi,bj)*drF(Nr)*recip_drF(ks)
     &                *hFacC(i,j,Nr,bi,bj)*recip_hFacC(i,j,ks,bi,bj)
              tmp_wS_M = tmp_wS_M
     &         + rA(i,j,bi,bj)*wVel(i,j,ks,bi,bj)*tmpVal
             ENDIF
            ENDDO
           ENDDO
           wS_Mean = tmp_wS_M
          ENDIF
#endif /* ALLOW_AIM */


C-- Atmos in Pot.Temp => conmpute energy conversion Temp -> PE 
C    = Omega*Theta*DeltaPI
          IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
           DO k=2,Nr
            ddPI=atm_cp*( (rC(K-1)/atm_po)**atm_kappa
     &                   -(rC( K )/atm_po)**atm_kappa )
            DO j=1,sNy
             DO i=1,sNx
              theta2PE = theta2PE
     &         - ddPI*rA(i,j,bi,bj)*wVel(i,j,k,bi,bj)
     &           *(theta(i,j,k,bi,bj)+theta(i,j,k-1,bi,bj))*0.5 _d 0
     &           *maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
             ENDDO
            ENDDO
           ENDDO
          ENDIF

#ifdef NONLIN_FRSURF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
        IF (select_rStar.NE.0) THEN 
C-- Compute Volume expansion term & do the integral
          vT_Mean = 0.
          vS_Mean = 0.
          vT_Heat = 0.
          DO k=1,Nr
           conv_th2Heat = atm_cp*((rC(k)/atm_po)**atm_kappa)
           DO j=1,sNy
            DO i=1,sNx
              tmpVol  = rA(i,j,bi,bj)*h0FacC(i,j,k,bi,bj)*drF(k)
              tmpVal  = rStarDhCDt(i,j,bi,bj)*theta(i,j,k,bi,bj)
              vT_Mean = vT_Mean + tmpVol*tmpVal
              vS_Mean = vS_Mean
     &          +tmpVol*rStarDhCDt(i,j,bi,bj)*salt(i,j,k,bi,bj)
C-  Atmos in Pot.Temp => convert Omega*Theta to heat flux :
              IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
               vT_Heat = vT_Heat + tmpVol*tmpVal*conv_th2Heat
              ENDIF
            ENDDO
           ENDDO
          ENDDO             
          wT_Mean = wT_Mean + vT_Mean
          wS_Mean = wS_Mean + vS_Mean
          wT_Heat = wT_Heat + vT_Heat
        ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* NONLIN_FRSURF */

C--    end bi,bj loop 
        ENDDO
       ENDDO

       _GLOBAL_SUM_R8(theArea,myThid) 
       _GLOBAL_SUM_R8(wT_Mean,myThid) 
       _GLOBAL_SUM_R8(wS_Mean,myThid) 
       IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
        _GLOBAL_SUM_R8(wT_Heat,myThid) 
        _GLOBAL_SUM_R8(theta2PE,myThid) 
       ENDIF
       IF (theArea.GT.0.) THEN
         wT_Mean = wT_Mean / theArea
         wS_Mean = wS_Mean / theArea
         wT_Heat = wT_Heat / theArea
         theta2PE = theta2PE / theArea
         wT_Heat  = wT_Heat  * rhoConst*recip_horiVertRatio
         theta2PE = theta2PE * rhoConst*recip_horiVertRatio
       ENDIF

C-    Print the Average value (monitor type output)

       CALL MON_SET_PREF('surfExpan',myThid)
       CALL MON_OUT_RL( '_theta', wT_Mean, mon_foot_mean ,myThid)
       CALL MON_OUT_RL( '_salt' , wS_Mean, mon_foot_mean ,myThid)
      IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
       CALL MON_OUT_RL( '_Heat' , wT_Heat, mon_foot_mean ,myThid)
       CALL MON_SET_PREF('En_Budget',myThid)
       CALL MON_OUT_RL('_T2PE',theta2PE, mon_foot_mean ,myThid)
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_surfcor.F,v 1.4 2003/12/24 00:29:17 jmc Exp $
2	C $Name: $
3
4	#include "MONITOR_OPTIONS.h"
5
6	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
7	CBOP
8	C !ROUTINE: MON_SURFCOR
9
10	C !INTERFACE:
11	SUBROUTINE MON_SURFCOR(
12	I myThid )
13
14	C !DESCRIPTION:
15	C Compute and write area-mean surface expansion term (also called
16	C ``surface correction'' with Linear FS).
17	C
18	C Diagnose mean surface expansion term
19	C \begin{equation}
20	C \mbox{with r coordinate} = (\mbox{w surf})(\mbox{Tracer})
21	C \end{equation}
22	C \begin{equation}
23	C \mbox{units} = (\mbox{W units})(\mbox{Tracer units})
24	C \ ; \ \mbox{+ = out}
25	C \end{equation}
26	C \begin{equation}
27	C \mbox{with r* coord}) = \frac{d\eta}{dt} \frac{dz}{H}
28	C (\mbox{Tracer})
29	C \end{equation}
30	C
31	C Atmosphere: convert surf.cor(Theta) to surface heating,
32	C \begin{equation}
33	C \mbox{units} = \frac{W}{m^2}, \mbox{+ = out}
34	C \end{equation}
35	C compute mean conversion term Temp -> PE , units= W/m2,
36	C + = decreasing heat content, increasing PE
37
38	C !USES:
39	IMPLICIT NONE
40	#include "SIZE.h"
41	#include "EEPARAMS.h"
42	#include "PARAMS.h"
43	#include "DYNVARS.h"
44	#include "SURFACE.h"
45	#include "GRID.h"
46	#include "MONITOR.h"
47
48	C !INPUT PARAMETERS:
49	INTEGER myThid
50	CEOP
51
52	C !LOCAL VARIABLES:
53	INTEGER i,j,k,ks,bi,bj
54	_RL theArea, wT_Mean, wS_Mean, tmp_wS_M, wT_Heat
55	_RL vT_Mean, vS_Mean, vT_Heat, theta2PE
56	_RL tmpVol, tmpVal, conv_th2Heat, ddPI
57
58	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
59	theArea = 0.
60	theta2PE = 0.
61	wT_Mean = 0.
62	wS_Mean = 0.
63	wT_Heat = 0.
64	vT_Mean = 0.
65	vS_Mean = 0.
66	vT_Heat = 0.
67	DO bj=myByLo(myThid),myByHi(myThid)
68	DO bi=myBxLo(myThid),myBxHi(myThid)
69	C-- Compute surface "expansion" term & do the integral
70	tmp_wS_M = wS_Mean
71	DO j=1,sNy
72	DO i=1,sNx
73	ks = ksurfC(i,j,bi,bj)
74	IF (ks.LE.Nr) THEN
75	theArea = theArea + rA(i,j,bi,bj)
76	tmpVal =
77	& rA(i,j,bi,bj)wVel(i,j,ks,bi,bj)theta(i,j,ks,bi,bj)
78	wT_Mean = wT_Mean + tmpVal
79	wS_Mean = wS_Mean
80	& + rA(i,j,bi,bj)wVel(i,j,ks,bi,bj) salt(i,j,ks,bi,bj)
81	C- Atmos in Pot.Temp => convert Omega*Theta to heat flux :
82	IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
83	wT_Heat = wT_Heat
84	& + tmpValatm_cp((rC(ks)/atm_po)**atm_kappa)
85	ENDIF
86	ENDIF
87	ENDDO
88	ENDDO
89	#ifdef ALLOW_AIM
90	IF ( useAIM ) THEN
91	DO j=1,sNy
92	DO i=1,sNx
93	ks = ksurfC(i,j,bi,bj)
94	IF (ks.LE.Nr) THEN
95	tmpVal = salt(i,j,ks,bi,bj)
96	& + salt(i,j,Nr,bi,bj)drF(Nr)recip_drF(ks)
97	& hFacC(i,j,Nr,bi,bj)recip_hFacC(i,j,ks,bi,bj)
98	tmp_wS_M = tmp_wS_M
99	& + rA(i,j,bi,bj)wVel(i,j,ks,bi,bj)tmpVal
100	ENDIF
101	ENDDO
102	ENDDO
103	wS_Mean = tmp_wS_M
104	ENDIF
105	#endif /* ALLOW_AIM */
106
107
108	C-- Atmos in Pot.Temp => conmpute energy conversion Temp -> PE
109	C = OmegaThetaDeltaPI
110	IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
111	DO k=2,Nr
112	ddPI=atm_cp( (rC(K-1)/atm_po)*atm_kappa
113	& -(rC( K )/atm_po)**atm_kappa )
114	DO j=1,sNy
115	DO i=1,sNx
116	theta2PE = theta2PE
117	& - ddPIrA(i,j,bi,bj)wVel(i,j,k,bi,bj)
118	& (theta(i,j,k,bi,bj)+theta(i,j,k-1,bi,bj))0.5 _d 0
119	& maskC(i,j,k-1,bi,bj)maskC(i,j,k,bi,bj)
120	ENDDO
121	ENDDO
122	ENDDO
123	ENDIF
124
125	#ifdef NONLIN_FRSURF
126	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
127	IF (select_rStar.NE.0) THEN
128	C-- Compute Volume expansion term & do the integral
129	vT_Mean = 0.
130	vS_Mean = 0.
131	vT_Heat = 0.
132	DO k=1,Nr
133	conv_th2Heat = atm_cp((rC(k)/atm_po)*atm_kappa)
134	DO j=1,sNy
135	DO i=1,sNx
136	tmpVol = rA(i,j,bi,bj)h0FacC(i,j,k,bi,bj)drF(k)
137	tmpVal = rStarDhCDt(i,j,bi,bj)*theta(i,j,k,bi,bj)
138	vT_Mean = vT_Mean + tmpVol*tmpVal
139	vS_Mean = vS_Mean
140	& +tmpVolrStarDhCDt(i,j,bi,bj)salt(i,j,k,bi,bj)
141	C- Atmos in Pot.Temp => convert Omega*Theta to heat flux :
142	IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
143	vT_Heat = vT_Heat + tmpVoltmpValconv_th2Heat
144	ENDIF
145	ENDDO
146	ENDDO
147	ENDDO
148	wT_Mean = wT_Mean + vT_Mean
149	wS_Mean = wS_Mean + vS_Mean
150	wT_Heat = wT_Heat + vT_Heat
151	ENDIF
152	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
153	#endif /* NONLIN_FRSURF */
154
155	C-- end bi,bj loop
156	ENDDO
157	ENDDO
158
159	_GLOBAL_SUM_R8(theArea,myThid)
160	_GLOBAL_SUM_R8(wT_Mean,myThid)
161	_GLOBAL_SUM_R8(wS_Mean,myThid)
162	IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
163	_GLOBAL_SUM_R8(wT_Heat,myThid)
164	_GLOBAL_SUM_R8(theta2PE,myThid)
165	ENDIF
166	IF (theArea.GT.0.) THEN
167	wT_Mean = wT_Mean / theArea
168	wS_Mean = wS_Mean / theArea
169	wT_Heat = wT_Heat / theArea
170	theta2PE = theta2PE / theArea
171	wT_Heat = wT_Heat * rhoConst*recip_horiVertRatio
172	theta2PE = theta2PE * rhoConst*recip_horiVertRatio
173	ENDIF
174
175	C- Print the Average value (monitor type output)
176
177	CALL MON_SET_PREF('surfExpan',myThid)
178	CALL MON_OUT_RL( '_theta', wT_Mean, mon_foot_mean ,myThid)
179	CALL MON_OUT_RL( '_salt' , wS_Mean, mon_foot_mean ,myThid)
180	IF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
181	CALL MON_OUT_RL( '_Heat' , wT_Heat, mon_foot_mean ,myThid)
182	CALL MON_SET_PREF('En_Budget',myThid)
183	CALL MON_OUT_RL('_T2PE',theta2PE, mon_foot_mean ,myThid)
184	ENDIF
185
186	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
187
188	RETURN
189	END