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

Contents of /MITgcm/pkg/my82/my82_ri_number.F

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

Revision 1.2 - (show annotations) (download)
Mon May 30 07:41:45 2005 UTC (19 years 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, checkpoint58r_post, checkpoint57i_post, checkpoint57y_post, checkpoint58n_post, checkpoint58x_post, checkpoint58t_post, checkpoint58h_post, checkpoint57y_pre, checkpoint58q_post, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, 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, checkpoint58p_post, checkpoint57j_post, checkpoint58b_post, checkpoint58m_post, checkpoint57l_post
Changes since 1.1: +2 -2 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 C $Header: /u/gcmpack/MITgcm/pkg/my82/my82_ri_number.F,v 1.1 2004/09/02 09:11:54 mlosch Exp $
2 C $Name: $
3
4 #include "MY82_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: MY82_RI_NUMBER
8
9 C !INTERFACE: ==========================================================
10 subroutine MY82_RI_NUMBER(
11 I bi, bj, K, iMin, iMax, jMin, jMax,
12 O RiNumber, buoyFreq, vertShear,
13 I myTime, myThid )
14
15 C !DESCRIPTION: \bv
16 C /==========================================================\
17 C | SUBROUTINE MY82_RI_NUMBER |
18 C | o Compute gradient Richardson number for Mellor and |
19 C | Yamada (1981) turbulence model |
20 C \==========================================================/
21 IMPLICIT NONE
22 c
23 c-------------------------------------------------------------------------
24
25 c \ev
26
27 C !USES: ===============================================================
28 #include "SIZE.h"
29 #include "EEPARAMS.h"
30 #include "PARAMS.h"
31 #include "DYNVARS.h"
32 #include "MY82.h"
33 #include "GRID.h"
34 #ifdef ALLOW_AUTODIFF_TAMC
35 #include "tamc.h"
36 #include "tamc_keys.h"
37 #else /* ALLOW_AUTODIFF_TAMC */
38 integer ikppkey
39 #endif /* ALLOW_AUTODIFF_TAMC */
40
41 C !INPUT PARAMETERS: ===================================================
42 C Routine arguments
43 C bi, bj - array indices on which to apply calculations
44 C iMin, iMax, jMin, jMax
45 C - array boundaries
46 C k - depth level
47 C myTime - Current time in simulation
48 C RiNumber - (output) Richardson number
49 C buoyFreq - (output) (neg.) buoyancy frequency -N^2
50 C vertShear - (output) vertical shear of velocity
51 INTEGER bi, bj, iMin, iMax, jMin, jMax, k
52 INTEGER myThid
53 _RL myTime
54 _RL RiNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55 _RL buoyFreq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56 _RL vertShear (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57
58 #ifdef ALLOW_MY82
59
60 C !LOCAL VARIABLES: ====================================================
61 C I,J,kUp - loop indices
62 C p0-125 - averaging coefficients
63 C tempu, tempv - temporary variables
64 C rhoK, rhoKm1 - Density below and above current interface
65 C epsilon - small number
66 INTEGER I,J,Km1
67 _RL p5 , p125
68 PARAMETER( p5=0.5D0, p125=0.125D0 )
69 _RL tempu, tempv
70 _RL epsilon
71 PARAMETER ( epsilon = 1.D-10 )
72 _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 _RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 #ifdef MY82_SMOOTH_RI
75 _RL RiTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 #endif /* MY82_SMOOTH_RI */
77 CEOP
78
79 Km1 = MAX(1,K-1)
80 C Not sure what is correct for pressure coordinates:
81 C RiFlux is always correct because it quadratic
82 C buoyFreq should also be correct in pressure coordinates:
83 C N^2=g^2drho/dp and K=1 at the bottom leads to the correct sign
84 C So the following is wrong:
85 CML IF ( buoyancyRelation .EQ. 'OCEANIC') THEN
86 CML kUp = MAX(1,K-1)
87 CML kDown = K
88 CML ELSEIF ( buoyancyRelation .EQ. 'OCEANIP') THEN
89 CML kUp = K
90 CML kDown = MAX(1,K-1)
91 CML ELSE
92 CML STOP 'MY82_RI_NUMBER: We should never reach this point'
93 CML ENDIF
94 C preparation: find density a kth and k-1st level
95 CALL FIND_RHO(
96 I bi, bj, iMin, iMax, jMin, jMax, Km1, K,
97 I theta, salt,
98 O rhoKm1,
99 I myThid )
100 CALL FIND_RHO(
101 I bi, bj, iMin, iMax, jMin, jMax, K, K,
102 I theta, salt,
103 O rhoK,
104 I myThid )
105
106 C first step: calculate flux Richardson number.
107 C calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.
108 DO J= jMin, jMax
109 DO I = iMin, iMax
110 tempu= .5*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
111 & - (uVel(I,J,K ,bi,bj)+uVel(I+1,J,K ,bi,bj)) )
112 & *recip_drC(K)
113 tempv= .5*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
114 & - (vVel(I,J,K ,bi,bj)+vVel(I,J+1,K ,bi,bj)) )
115 & *recip_drC(K)
116 vertShear(I,J) = tempu*tempu+tempv*tempv
117
118 C
119 C calculate g*(drho/dz)/rho0= -N^2 .
120 C
121 buoyFreq(I,J) = gravity*recip_rhoConst*horiVertRatio *
122 & (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)
123 C
124 C calculates gradient Richardson number, bounded by
125 C a very large negative value.
126 C
127 RiNumber(I,J) = -buoyFreq(I,J)/max(vertShear(I,J),epsilon)
128
129 ENDDO
130 ENDDO
131
132 #ifdef MY82_SMOOTH_RI
133 C average Richardson number horizontally
134 DO J=jMin,jMax
135 DO I=iMin,iMax
136 RiTmp(I,J) = RiNumber(I,J)
137 ENDDO
138 ENDDO
139 DO J=1-Oly+1,sNy+Oly-1
140 DO I=1-Olx+1,sNx+Olx-1
141 RiNumber(I,J) = p5*RiTmp(I,J)
142 & + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J)
143 & + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1)
144 ENDDO
145 ENDDO
146 #endif /* MY82_SMOOTH_RI */
147
148 #endif /* ALLOW_MY82 */
149
150 RETURN
151 END
152
  ViewVC Help
Powered by ViewVC 1.1.22