/[MITgcm]/MITgcm/pkg/my82/my82_calc.F
ViewVC logotype

Annotation of /MITgcm/pkg/my82/my82_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.4 - (hide annotations) (download)
Mon May 30 07:41:45 2005 UTC (18 years, 11 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint57t_post, checkpoint57o_post, checkpoint58e_post, checkpoint57v_post, checkpoint58u_post, checkpoint58w_post, checkpoint57m_post, checkpoint57s_post, checkpoint57k_post, checkpoint60, checkpoint61, checkpoint58r_post, checkpoint57i_post, checkpoint57y_post, checkpoint58n_post, checkpoint58x_post, checkpoint58t_post, checkpoint58h_post, checkpoint57y_pre, checkpoint58q_post, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint57r_post, checkpoint59, checkpoint58, checkpoint58f_post, checkpoint57x_post, checkpoint57n_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint58y_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint58p_post, checkpoint61a, checkpoint57j_post, checkpoint58b_post, checkpoint58m_post, checkpoint57l_post
Changes since 1.3: +10 -10 lines 
added a few "_d 0" and D0 (in PARAMETER statements) in a desperate
effort to make vermix.my82 pass on a Sun, unfortunately no success;
but for cleaner code I check it in anyway.
1 mlosch 1.4 C $Header: /u/gcmpack/MITgcm/pkg/my82/my82_calc.F,v 1.3 2005/04/06 18:43:11 jmc 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     #include "DYNVARS.h"
35     #include "MY82.h"
36     #include "FFIELDS.h"
37     #include "GRID.h"
38     #ifdef ALLOW_AUTODIFF_TAMC
39     #include "tamc.h"
40     #include "tamc_keys.h"
41     #else /* ALLOW_AUTODIFF_TAMC */
42     integer ikppkey
43     #endif /* ALLOW_AUTODIFF_TAMC */
44    
45     C !INPUT PARAMETERS: ===================================================
46     c Routine arguments
47     c bi, bj - array indices on which to apply calculations
48     c myTime - Current time in simulation
49    
50     INTEGER bi, bj
51     INTEGER myThid
52     _RL myTime
53    
54     #ifdef ALLOW_MY82
55    
56     C !LOCAL VARIABLES: ====================================================
57     c Local constants
58     C imin, imax, jmin, jmax - array computation indices
59     C RiNumber - Richardson Number = -GH/GM
60     C GH - buoyancy freqency after call to Ri_number, later
61     C GH of M. Satoh, Eq. (11.3.45)
62     C GM - vertical shear of velocity after call Ri_number,
63     C later GM of M. Satoh, Eq. (11.3.44)
64     INTEGER I, J, K
65     INTEGER iMin ,iMax ,jMin ,jMax
66     _RL RiFlux, RiTmp
67     _RL SHtmp, bTmp, tkesquare, tkel
68     _RL RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69     _RL GH(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70     _RL GM(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RL SH(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
72     _RL SM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73     _RL tke(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
74     CEOP
75     iMin = 2-OLx
76     iMax = sNx+OLx-1
77     jMin = 2-OLy
78     jMax = sNy+OLy-1
79    
80     C Initialize local fields
81     DO J=1-Oly,sNy+Oly
82     DO I=1-Olx,sNx+Olx
83 mlosch 1.4 GH(I,J) = 0. _d 0
84     GM(I,J) = 0. _d 0
85 mlosch 1.1 ENDDO
86     ENDDO
87     DO K = 1, Nr
88     DO J=1-Oly,sNy+Oly
89     DO I=1-Olx,sNx+Olx
90 mlosch 1.4 SH(I,J,K) = 0. _d 0
91     SM(I,J,K) = 0. _d 0
92     tke(I,J,K) = 0. _d 0
93 mlosch 1.1 ENDDO
94     ENDDO
95     ENDDO
96     C first k-loop
97     C compute turbulent kinetic energy from richardson number
98     DO K = 2, Nr
99     CALL MY82_RI_NUMBER(
100     I bi, bj, K, iMin, iMax, jMin, jMax,
101     O RiNumber, GH, GM,
102     I myTime, myThid )
103     DO J=jMin,jMax
104     DO I=iMin,iMax
105     RiTmp = MIN(RiNumber(I,J),RiMax)
106     btmp = beta1+beta4*RiTmp
107     C M. Satoh, Atmospheric Circulation Dynamics and General
108     C Circulation models, Springer, 2004: Eq. (11.3.60)
109     RiFlux = ( btmp - SQRT(btmp*btmp - 4.*beta2*beta3*RiTmp) )
110     & /(2.*beta2)
111     C M. Satoh: Eq. (11.3.58)
112     SHtmp = (alpha1-alpha2*RiFlux)/(1.-RiFlux)
113     SH(I,J,K) = SHtmp
114     SM(I,J,K) = SHtmp*(beta1-beta2*RiFlux)/(beta3-beta4*RiFlux)
115     C M. Satoh: Eq. (11.3.53/55)
116     tkesquare = MAX( 0. _d 0,
117     & b1*(SH(I,J,K)*GH(I,J) + SM(I,J,K)*GM(I,J)) )
118     tke(I,J,K) = sqrt(tkesquare)
119     CML if ( k.eq.2 .and. i.ge.1 .and. i.le.sNx .and. j.eq.1)
120     CML & print '(A,3I3,8E13.5)', 'ml-my82', I,J,K, RiNumber(I,J),
121     CML & RiTmp,RiFlux,
122     CML & SH(I,J,K), SM(I,J,K), GH(I,J), GM(I,J), tke(I,J,K)
123     ENDDO
124     ENDDO
125     C end of first k-loop
126     ENDDO
127    
128     C re-initilialize GM and GH for abuse, they no longer have
129     C the meaning of shear and negative buoyancy frequency
130     DO J=jMin,jMax
131     DO I=iMin,iMax
132 mlosch 1.4 GH(I,J) = 0. _d 0
133     GM(I,J) = 0. _d 0
134 mlosch 1.1 ENDDO
135     ENDDO
136     C Find boundary length scale from energy weighted mean.
137     C This is something like "center of mass" of the vertical
138     C tke-distribution
139     C begin second k-loop
140     DO K = 2, Nr
141     DO J=jMin,jMax
142     DO I=iMin,iMax
143     GM(I,J) = GM(I,J) + tke(I,J,K)*rF(K)
144     GH(I,J) = GH(I,J) + tke(I,J,K)
145     ENDDO
146     ENDDO
147     C end of second k-loop
148     END DO
149     C compute boundary length scale MYhbl
150     DO J=jMin,jMax
151     DO I=iMin,iMax
152 mlosch 1.4 IF ( GH(I,J) .EQ. 0. _d 0 ) THEN
153     MYhbl(I,J,bi,bj) = 0. _d 0
154 mlosch 1.1 ELSE
155     MYhbl(I,J,bi,bj) = -GM(I,J)/GH(I,J)*MYhblScale
156     ENDIF
157     ENDDO
158     ENDDO
159     C begin third k-loop
160     DO K = 1, Nr
161     C integrate tke to find integral length scale
162     DO J=jMin,jMax
163     DO I=iMin,iMax
164     tkel = MYhbl(I,J,bi,bj)*tke(I,J,K)
165     C M. Satoh: Eq. (11.3.43)
166     MYviscAr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)*tkel*SM(I,J,K)
167     MYdiffKr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)*tkel*SH(I,J,K)
168     C Set a minium (= background) value
169     MYviscAr(I,J,K,bi,bj) = MAX(MYviscAr(I,J,K,bi,bj),viscAr)
170 jmc 1.2 MYdiffKr(I,J,K,bi,bj) = MAX(MYdiffKr(I,J,K,bi,bj),
171     & diffKrNrT(k) )
172 mlosch 1.1 C Set a maximum and mask land point
173     MYviscAr(I,J,K,bi,bj) = MIN(MYviscAr(I,J,K,bi,bj),MYviscMax)
174     & * maskC(I,J,K,bi,bj)
175     MYdiffKr(I,J,K,bi,bj) = MIN(MYdiffKr(I,J,K,bi,bj),MYdiffMax)
176     & * maskC(I,J,K,bi,bj)
177     ENDDO
178     ENDDO
179     C end third k-loop
180     ENDDO
181    
182     #endif /* ALLOW_MY82 */
183    
184     RETURN
185     END
186    
  ViewVC Help
Powered by ViewVC 1.1.22