pkg/monitor/mon_surfcor.F

C $Header: /u/gcmpack/MITgcm/pkg/monitor/mon_surfcor.F,v 1.1 2002/11/22 02:33:24 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

      SUBROUTINE MON_SURFCOR(
     I                       myThid )
C     |==========================================================
C     | SUBROUTINE MON_SURFCOR                                   
C     | o Compute and write area-mean surface expansion term    
C     |   (also called "surface correction" with Linear FS)
C     |==========================================================
C     | Diagnose mean surface expansion term
C     |             with r coordinate  = w_surf * Tracer
C     |                 units = W_units * Tracer units ; + = out 
C     |             with r* coord.:  = d.eta/dt * dz/H * Tracer
C     | Atmosphere :                                               
C     |   convert surf.cor(Theta) to surface heating, 
C     |                 units= W/m2, + = out
C     |   compute mean conversion term Temp -> PE , units= W/m2, 
C     |               + = decreasing heat content, increasing PE
C     |==========================================================
      IMPLICIT NONE

C     === Global data ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "SURFACE.h"
#include "GRID.h"
#include "MONITOR.h"

C     === Routine arguments ===
      INTEGER myThid

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