pkg/my82/my82_calc.F

C $Header: /u/gcmpack/MITgcm/pkg/my82/my82_calc.F,v 1.4 2005/05/30 07:41:45 mlosch Exp $
C $Name:  $

#include "MY82_OPTIONS.h"

CBOP
C !ROUTINE: MY82_CALC

C !INTERFACE: ======================================================
      subroutine MY82_CALC(
     I     bi, bj, myTime, myThid )

C !DESCRIPTION: \bv
C     /==========================================================\
C     | SUBROUTINE MY82_CALC                                     |
C     | o Compute all MY82 fields defined in MY82.h              |
C     |==========================================================|
C     | This subroutine is based on SPEM code                    |
C     \==========================================================/
      IMPLICIT NONE
C
C--------------------------------------------------------------------

C global parameters updated by pp_calc
C     PPviscAz   - PP eddy viscosity coefficient              (m^2/s)
C     PPdiffKzT  - PP diffusion coefficient for temperature   (m^2/s)
C
C \ev

C !USES: ============================================================
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
c#include "DYNVARS.h"
#include "MY82.h"
c#include "FFIELDS.h"
#include "GRID.h"

C !INPUT PARAMETERS: ===================================================
c Routine arguments
c     bi, bj - array indices on which to apply calculations
c     myTime - Current time in simulation

      INTEGER bi, bj
      _RL     myTime
      INTEGER myThid

#ifdef ALLOW_MY82

C !LOCAL VARIABLES: ====================================================
c Local constants
C     imin, imax, jmin, jmax  - array computation indices
C     RiNumber - Richardson Number = -GH/GM
C     GH       - buoyancy freqency after call to Ri_number, later
C                GH of M. Satoh, Eq. (11.3.45)
C     GM       - vertical shear of velocity after call Ri_number,
C                later GM of M. Satoh, Eq. (11.3.44)
      INTEGER I, J, K
      INTEGER   iMin ,iMax ,jMin ,jMax
      _RL     RiFlux, RiTmp
      _RL     SHtmp, bTmp, tkesquare, tkel
      _RL     RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     GH(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     GM(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     SH(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     SM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     tke(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
CEOP
      iMin = 2-OLx
      iMax = sNx+OLx-1
      jMin = 2-OLy
      jMax = sNy+OLy-1

C     Initialize local fields
      DO J=1-Oly,sNy+Oly
       DO I=1-Olx,sNx+Olx
        GH(I,J) = 0. _d 0
        GM(I,J) = 0. _d 0
       ENDDO
      ENDDO
      DO K = 1, Nr
       DO J=1-Oly,sNy+Oly
        DO I=1-Olx,sNx+Olx
         SH(I,J,K)  = 0. _d 0
         SM(I,J,K)  = 0. _d 0
         tke(I,J,K) = 0. _d 0
        ENDDO
       ENDDO
      ENDDO
C     first k-loop
C     compute turbulent kinetic energy from richardson number
      DO K = 2, Nr
       CALL MY82_RI_NUMBER(
     I      bi, bj, K, iMin, iMax, jMin, jMax,
     O      RiNumber, GH, GM,
     I      myTime, myThid )
       DO J=jMin,jMax
        DO I=iMin,iMax
         RiTmp  = MIN(RiNumber(I,J),RiMax)
         btmp   = beta1+beta4*RiTmp
C     M. Satoh, Atmospheric Circulation Dynamics and General
C     Circulation models, Springer, 2004: Eq. (11.3.60)
         RiFlux =( btmp - SQRT(btmp*btmp - 4. _d 0 *beta2*beta3*RiTmp) )
     &        /(2. _d 0*beta2)
C     M. Satoh: Eq. (11.3.58)
         SHtmp       = (alpha1-alpha2*RiFlux)/(1. _d 0-RiFlux)
         SH(I,J,K)   = SHtmp
         SM(I,J,K)   = SHtmp*(beta1-beta2*RiFlux)/(beta3-beta4*RiFlux)
C     M. Satoh: Eq. (11.3.53/55)
         tkesquare = MAX( 0. _d 0,
     &        b1*(SH(I,J,K)*GH(I,J) + SM(I,J,K)*GM(I,J)) )
         tke(I,J,K) = sqrt(tkesquare)
CML         if ( k.eq.2 .and. i.ge.1 .and. i.le.sNx .and. j.eq.1)
CML     &   print '(A,3I3,8E13.5)', 'ml-my82', I,J,K, RiNumber(I,J),
CML     &        RiTmp,RiFlux,
CML     &      SH(I,J,K), SM(I,J,K), GH(I,J), GM(I,J), tke(I,J,K)
        ENDDO
       ENDDO
C     end of first k-loop
      ENDDO

C     re-initilialize GM and GH for abuse, they no longer have
C     the meaning of shear and negative  buoyancy frequency
      DO J=jMin,jMax
       DO I=iMin,iMax
        GH(I,J) = 0. _d 0
        GM(I,J) = 0. _d 0
       ENDDO
      ENDDO
C     Find boundary length scale from energy weighted mean.
C     This is something like "center of mass" of the vertical
C     tke-distribution
C     begin second k-loop
      DO K = 2, Nr
       DO J=jMin,jMax
        DO I=iMin,iMax
         GM(I,J) = GM(I,J) + tke(I,J,K)*rF(K)
         GH(I,J) = GH(I,J) + tke(I,J,K)
        ENDDO
       ENDDO
C     end of second k-loop
      END DO
C     compute boundary length scale MYhbl
      DO J=jMin,jMax
       DO I=iMin,iMax
        IF ( GH(I,J) .EQ. 0. _d 0 ) THEN
         MYhbl(I,J,bi,bj) = 0. _d 0
        ELSE
         MYhbl(I,J,bi,bj) = -GM(I,J)/GH(I,J)*MYhblScale
        ENDIF
       ENDDO
      ENDDO
C     begin third k-loop
      DO K = 1, Nr
C     integrate tke to find integral length scale
       DO J=jMin,jMax
        DO I=iMin,iMax
         tkel = MYhbl(I,J,bi,bj)*tke(I,J,K)
C     M. Satoh: Eq. (11.3.43)
         MYviscAr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)*tkel*SM(I,J,K)
         MYdiffKr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)*tkel*SH(I,J,K)
C     Set a minium (= background) value
         MYviscAr(I,J,K,bi,bj) = MAX(MYviscAr(I,J,K,bi,bj),viscAr)
         MYdiffKr(I,J,K,bi,bj) = MAX(MYdiffKr(I,J,K,bi,bj),
     &                               diffKrNrT(k) )
C     Set a maximum and mask land point
         MYviscAr(I,J,K,bi,bj) = MIN(MYviscAr(I,J,K,bi,bj),MYviscMax)
     &        * maskC(I,J,K,bi,bj)
         MYdiffKr(I,J,K,bi,bj) = MIN(MYdiffKr(I,J,K,bi,bj),MYdiffMax)
     &        * maskC(I,J,K,bi,bj)
        ENDDO
       ENDDO
C     end third k-loop
      ENDDO

#endif /* ALLOW_MY82 */

      RETURN
      END
1	jmc	1.5	C $Header: /u/gcmpack/MITgcm/pkg/my82/my82_calc.F,v 1.4 2005/05/30 07:41:45 mlosch Exp $
2	mlosch	1.1	C $Name: $
3
4			#include "MY82_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: MY82_CALC
8
9			C !INTERFACE: ======================================================
10			subroutine MY82_CALC(
11			I bi, bj, myTime, myThid )
12
13			C !DESCRIPTION: \bv
14			C /==========================================================\
15			C \| SUBROUTINE MY82_CALC \|
16			C \| o Compute all MY82 fields defined in MY82.h \|
17			C \|==========================================================\|
18			C \| This subroutine is based on SPEM code \|
19			C \==========================================================/
20			IMPLICIT NONE
21			C
22			C--------------------------------------------------------------------
23
24			C global parameters updated by pp_calc
25			C PPviscAz - PP eddy viscosity coefficient (m^2/s)
26			C PPdiffKzT - PP diffusion coefficient for temperature (m^2/s)
27			C
28			C \ev
29
30			C !USES: ============================================================
31			#include "SIZE.h"
32			#include "EEPARAMS.h"
33			#include "PARAMS.h"
34	jmc	1.5	c#include "DYNVARS.h"
35	mlosch	1.1	#include "MY82.h"
36	jmc	1.5	c#include "FFIELDS.h"
37	mlosch	1.1	#include "GRID.h"
38
39			C !INPUT PARAMETERS: ===================================================
40			c Routine arguments
41			c bi, bj - array indices on which to apply calculations
42			c myTime - Current time in simulation
43
44			INTEGER bi, bj
45	jmc	1.5	_RL myTime
46	mlosch	1.1	INTEGER myThid
47
48			#ifdef ALLOW_MY82
49
50			C !LOCAL VARIABLES: ====================================================
51			c Local constants
52			C imin, imax, jmin, jmax - array computation indices
53			C RiNumber - Richardson Number = -GH/GM
54			C GH - buoyancy freqency after call to Ri_number, later
55			C GH of M. Satoh, Eq. (11.3.45)
56			C GM - vertical shear of velocity after call Ri_number,
57			C later GM of M. Satoh, Eq. (11.3.44)
58			INTEGER I, J, K
59			INTEGER iMin ,iMax ,jMin ,jMax
60			_RL RiFlux, RiTmp
61			_RL SHtmp, bTmp, tkesquare, tkel
62			_RL RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63			_RL GH(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64			_RL GM(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65			_RL SH(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
66			_RL SM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
67			_RL tke(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
68			CEOP
69			iMin = 2-OLx
70			iMax = sNx+OLx-1
71			jMin = 2-OLy
72			jMax = sNy+OLy-1
73
74			C Initialize local fields
75			DO J=1-Oly,sNy+Oly
76			DO I=1-Olx,sNx+Olx
77	jmc	1.5	GH(I,J) = 0. _d 0
78			GM(I,J) = 0. _d 0
79	mlosch	1.1	ENDDO
80			ENDDO
81			DO K = 1, Nr
82			DO J=1-Oly,sNy+Oly
83			DO I=1-Olx,sNx+Olx
84	mlosch	1.4	SH(I,J,K) = 0. _d 0
85			SM(I,J,K) = 0. _d 0
86			tke(I,J,K) = 0. _d 0
87	mlosch	1.1	ENDDO
88	jmc	1.5	ENDDO
89	mlosch	1.1	ENDDO
90			C first k-loop
91			C compute turbulent kinetic energy from richardson number
92			DO K = 2, Nr
93			CALL MY82_RI_NUMBER(
94	jmc	1.5	I bi, bj, K, iMin, iMax, jMin, jMax,
95	mlosch	1.1	O RiNumber, GH, GM,
96			I myTime, myThid )
97			DO J=jMin,jMax
98			DO I=iMin,iMax
99			RiTmp = MIN(RiNumber(I,J),RiMax)
100			btmp = beta1+beta4*RiTmp
101			C M. Satoh, Atmospheric Circulation Dynamics and General
102			C Circulation models, Springer, 2004: Eq. (11.3.60)
103	jmc	1.5	RiFlux =( btmp - SQRT(btmpbtmp - 4. _d 0 beta2beta3RiTmp) )
104			& /(2. _d 0*beta2)
105	mlosch	1.1	C M. Satoh: Eq. (11.3.58)
106	jmc	1.5	SHtmp = (alpha1-alpha2*RiFlux)/(1. _d 0-RiFlux)
107	mlosch	1.1	SH(I,J,K) = SHtmp
108			SM(I,J,K) = SHtmp(beta1-beta2RiFlux)/(beta3-beta4*RiFlux)
109			C M. Satoh: Eq. (11.3.53/55)
110			tkesquare = MAX( 0. _d 0,
111			& b1(SH(I,J,K)GH(I,J) + SM(I,J,K)*GM(I,J)) )
112			tke(I,J,K) = sqrt(tkesquare)
113			CML if ( k.eq.2 .and. i.ge.1 .and. i.le.sNx .and. j.eq.1)
114			CML & print '(A,3I3,8E13.5)', 'ml-my82', I,J,K, RiNumber(I,J),
115			CML & RiTmp,RiFlux,
116			CML & SH(I,J,K), SM(I,J,K), GH(I,J), GM(I,J), tke(I,J,K)
117			ENDDO
118			ENDDO
119			C end of first k-loop
120	jmc	1.5	ENDDO
121	mlosch	1.1
122	jmc	1.5	C re-initilialize GM and GH for abuse, they no longer have
123	mlosch	1.1	C the meaning of shear and negative buoyancy frequency
124			DO J=jMin,jMax
125			DO I=iMin,iMax
126	mlosch	1.4	GH(I,J) = 0. _d 0
127			GM(I,J) = 0. _d 0
128	mlosch	1.1	ENDDO
129			ENDDO
130			C Find boundary length scale from energy weighted mean.
131			C This is something like "center of mass" of the vertical
132			C tke-distribution
133			C begin second k-loop
134			DO K = 2, Nr
135			DO J=jMin,jMax
136			DO I=iMin,iMax
137			GM(I,J) = GM(I,J) + tke(I,J,K)*rF(K)
138	jmc	1.5	GH(I,J) = GH(I,J) + tke(I,J,K)
139	mlosch	1.1	ENDDO
140			ENDDO
141			C end of second k-loop
142			END DO
143			C compute boundary length scale MYhbl
144			DO J=jMin,jMax
145			DO I=iMin,iMax
146	mlosch	1.4	IF ( GH(I,J) .EQ. 0. _d 0 ) THEN
147			MYhbl(I,J,bi,bj) = 0. _d 0
148	mlosch	1.1	ELSE
149			MYhbl(I,J,bi,bj) = -GM(I,J)/GH(I,J)*MYhblScale
150			ENDIF
151			ENDDO
152	jmc	1.5	ENDDO
153	mlosch	1.1	C begin third k-loop
154			DO K = 1, Nr
155			C integrate tke to find integral length scale
156			DO J=jMin,jMax
157			DO I=iMin,iMax
158			tkel = MYhbl(I,J,bi,bj)*tke(I,J,K)
159			C M. Satoh: Eq. (11.3.43)
160			MYviscAr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)tkelSM(I,J,K)
161			MYdiffKr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)tkelSH(I,J,K)
162			C Set a minium (= background) value
163			MYviscAr(I,J,K,bi,bj) = MAX(MYviscAr(I,J,K,bi,bj),viscAr)
164	jmc	1.2	MYdiffKr(I,J,K,bi,bj) = MAX(MYdiffKr(I,J,K,bi,bj),
165			& diffKrNrT(k) )
166	mlosch	1.1	C Set a maximum and mask land point
167			MYviscAr(I,J,K,bi,bj) = MIN(MYviscAr(I,J,K,bi,bj),MYviscMax)
168			& * maskC(I,J,K,bi,bj)
169			MYdiffKr(I,J,K,bi,bj) = MIN(MYdiffKr(I,J,K,bi,bj),MYdiffMax)
170			& * maskC(I,J,K,bi,bj)
171			ENDDO
172	jmc	1.5	ENDDO
173	mlosch	1.1	C end third k-loop
174	jmc	1.5	ENDDO
175	mlosch	1.1
176			#endif /* ALLOW_MY82 */
177
178			RETURN
179			END