1 | MODULE ldfc1d_c2d |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE ldfc1d_c2d *** |
---|
4 | !! Ocean physics: profile and horizontal shape of lateral eddy coefficients |
---|
5 | !!===================================================================== |
---|
6 | !! History : 3.7 ! 2013-12 (G. Madec) restructuration/simplification of aht/aeiv specification, |
---|
7 | !! ! add velocity dependent coefficient and optional read in file |
---|
8 | !!---------------------------------------------------------------------- |
---|
9 | |
---|
10 | !!---------------------------------------------------------------------- |
---|
11 | !! ldf_c1d : ah reduced by 1/4 on the vertical (tanh profile, inflection at 300m) |
---|
12 | !! ldf_c2d : ah = F(e1,e2) (laplacian or = F(e1^3,e2^3) (bilaplacian) |
---|
13 | !!---------------------------------------------------------------------- |
---|
14 | USE oce ! ocean dynamics and tracers |
---|
15 | USE dom_oce ! ocean space and time domain |
---|
16 | USE phycst ! physical constants |
---|
17 | ! |
---|
18 | USE in_out_manager ! I/O manager |
---|
19 | USE lib_mpp ! distribued memory computing library |
---|
20 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
21 | |
---|
22 | IMPLICIT NONE |
---|
23 | PRIVATE |
---|
24 | |
---|
25 | PUBLIC ldf_c1d ! called by ldftra and ldfdyn modules |
---|
26 | PUBLIC ldf_c2d ! called by ldftra and ldfdyn modules |
---|
27 | |
---|
28 | |
---|
29 | !! * Substitutions |
---|
30 | # include "vectopt_loop_substitute.h90" |
---|
31 | !!---------------------------------------------------------------------- |
---|
32 | !! NEMO/OPA 3.7 , NEMO Consortium (2015) |
---|
33 | !! $Id: $ |
---|
34 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
---|
35 | !!---------------------------------------------------------------------- |
---|
36 | CONTAINS |
---|
37 | |
---|
38 | SUBROUTINE ldf_c1d( cd_type, prat, pahs1, pahs2, pah1, pah2 ) |
---|
39 | !!---------------------------------------------------------------------- |
---|
40 | !! *** ROUTINE ldf_c1d *** |
---|
41 | !! |
---|
42 | !! ** Purpose : 1D eddy diffusivity/viscosity coefficients |
---|
43 | !! |
---|
44 | !! ** Method : 1D eddy diffusivity coefficients F( depth ) |
---|
45 | !! Reduction by prat from surface to bottom |
---|
46 | !! hyperbolic tangent profile with inflection point |
---|
47 | !! at zh=500m and a width of zw=200m |
---|
48 | !! |
---|
49 | !! cd_type = TRA pah1, pah2 defined at U- and V-points |
---|
50 | !! DYN pah1, pah2 defined at T- and F-points |
---|
51 | !!---------------------------------------------------------------------- |
---|
52 | CHARACTER(len=2) , INTENT(in ) :: cd_type ! DYNamique or TRAcers |
---|
53 | REAL(wp) , INTENT(in ) :: prat ! ratio surface/deep values [-] |
---|
54 | REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pahs1, pahs2 ! surface value of eddy coefficient [m2/s] |
---|
55 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pah1 , pah2 ! eddy coefficient [m2/s] |
---|
56 | ! |
---|
57 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
58 | REAL(wp) :: zh, zc, zdep1 ! local scalars |
---|
59 | REAL(wp) :: zw , zdep2 ! - - |
---|
60 | !!---------------------------------------------------------------------- |
---|
61 | |
---|
62 | IF(lwp) THEN |
---|
63 | WRITE(numout,*) |
---|
64 | WRITE(numout,*) ' ldf_c1d : set a given profile to eddy diffusivity/viscosity coefficients' |
---|
65 | WRITE(numout,*) ' ~~~~~~~' |
---|
66 | ENDIF |
---|
67 | |
---|
68 | ! initialization of the profile |
---|
69 | zh = 500._wp ! depth of the inflection point [m] |
---|
70 | zw = 1._wp / 200._wp ! width^-1 - - - [1/m] |
---|
71 | ! ! associated coefficient [-] |
---|
72 | zc = ( 1._wp - prat ) / ( 1._wp + TANH( zh * zw) ) |
---|
73 | ! |
---|
74 | ! |
---|
75 | SELECT CASE( cd_type ) ! point of calculation |
---|
76 | ! |
---|
77 | CASE( 'DYN' ) ! T- and F-points |
---|
78 | DO jk = 1, jpk ! pah1 at T-point |
---|
79 | pah1(:,:,jk) = pahs1(:,:) * ( prat + zc * ( 1._wp + TANH( - ( gdept_n(:,:,jk) - zh ) * zw) ) ) * tmask(:,:,jk) |
---|
80 | END DO |
---|
81 | DO jk = 1, jpk ! pah2 at F-point (zdep2 is an approximation in zps-coord.) |
---|
82 | DO jj = 1, jpjm1 |
---|
83 | DO ji = 1, fs_jpim1 |
---|
84 | zdep2 = ( gdept_n(ji,jj+1,jk) + gdept_n(ji+1,jj+1,jk) & |
---|
85 | & + gdept_n(ji,jj ,jk) + gdept_n(ji+1,jj ,jk) ) * 0.25_wp |
---|
86 | pah2(ji,jj,jk) = pahs2(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) * fmask(ji,jj,jk) |
---|
87 | END DO |
---|
88 | END DO |
---|
89 | END DO |
---|
90 | CALL lbc_lnk( pah2, 'F', 1. ) ! Lateral boundary conditions |
---|
91 | ! |
---|
92 | CASE( 'TRA' ) ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.) |
---|
93 | DO jk = 1, jpk |
---|
94 | DO jj = 1, jpjm1 |
---|
95 | DO ji = 1, fs_jpim1 |
---|
96 | zdep1 = ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) * 0.5_wp |
---|
97 | zdep2 = ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) * 0.5_wp |
---|
98 | pah1(ji,jj,jk) = pahs1(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) * umask(ji,jj,jk) |
---|
99 | pah2(ji,jj,jk) = pahs2(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) * vmask(ji,jj,jk) |
---|
100 | END DO |
---|
101 | END DO |
---|
102 | END DO |
---|
103 | CALL lbc_lnk( pah1, 'U', 1. ) ! Lateral boundary conditions |
---|
104 | CALL lbc_lnk( pah2, 'V', 1. ) |
---|
105 | ! |
---|
106 | END SELECT |
---|
107 | ! |
---|
108 | END SUBROUTINE ldf_c1d |
---|
109 | |
---|
110 | |
---|
111 | SUBROUTINE ldf_c2d( cd_type, cd_op, pah0, pah1, pah2 ) |
---|
112 | !!---------------------------------------------------------------------- |
---|
113 | !! *** ROUTINE ldf_c2d *** |
---|
114 | !! |
---|
115 | !! ** Purpose : 2D eddy diffusivity/viscosity coefficients |
---|
116 | !! |
---|
117 | !! ** Method : 2D eddy diffusivity coefficients F( e1 , e2 ) |
---|
118 | !! laplacian operator : ah proportional to the scale factor [m2/s] |
---|
119 | !! bilaplacian operator : ah proportional to the (scale factor)^3 [m4/s] |
---|
120 | !! In both cases, pah0 is the maximum value reached by the coefficient |
---|
121 | !! at the Equator in case of e1=ra*rad= ~111km, not over the whole domain. |
---|
122 | !! |
---|
123 | !! cd_type = TRA pah1, pah2 defined at U- and V-points |
---|
124 | !! DYN pah1, pah2 defined at T- and F-points |
---|
125 | !!---------------------------------------------------------------------- |
---|
126 | CHARACTER(len=3) , INTENT(in ) :: cd_type ! DYNamique or TRAcers |
---|
127 | CHARACTER(len=3) , INTENT(in ) :: cd_op ! operator: LAPlacian BiLaPlacian |
---|
128 | REAL(wp) , INTENT(in ) :: pah0 ! eddy coefficient [m2/s or m4/s] |
---|
129 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: pah1, pah2 ! eddy coefficient [m2/s or m4/s] |
---|
130 | ! |
---|
131 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
132 | REAL(wp) :: za00, zd_max, zemax1, zemax2 ! local scalar |
---|
133 | !!---------------------------------------------------------------------- |
---|
134 | ! |
---|
135 | zd_max = ra * rad ! = 1 degree at the equator in meters |
---|
136 | ! |
---|
137 | IF(lwp) THEN |
---|
138 | WRITE(numout,*) |
---|
139 | WRITE(numout,*) ' ldf_c2d : aht = rn_aht0 * max(e1,e2)/e_equator ( laplacian) ' |
---|
140 | WRITE(numout,*) ' ~~~~~~~ or = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)' |
---|
141 | WRITE(numout,*) |
---|
142 | ENDIF |
---|
143 | ! |
---|
144 | SELECT CASE( cd_type ) !== surface values ==! (depending on DYN/TRA) |
---|
145 | ! |
---|
146 | CASE( 'DYN' ) ! T- and F-points |
---|
147 | IF( cd_op == 'LAP' ) THEN ! laplacian operator |
---|
148 | IF(lwp) WRITE(numout,*) ' momentum laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' |
---|
149 | za00 = pah0 / zd_max |
---|
150 | DO jj = 1, jpj |
---|
151 | DO ji = 1, jpi |
---|
152 | zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) |
---|
153 | zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) |
---|
154 | pah1(ji,jj,1) = za00 * zemax1 |
---|
155 | pah2(ji,jj,1) = za00 * zemax2 |
---|
156 | END DO |
---|
157 | END DO |
---|
158 | ELSEIF( cd_op == 'BLP' ) THEN ! bilaplacian operator |
---|
159 | IF(lwp) WRITE(numout,*) ' momentum bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' |
---|
160 | za00 = pah0 / ( zd_max * zd_max * zd_max ) |
---|
161 | DO jj = 1, jpj |
---|
162 | DO ji = 1, jpi |
---|
163 | zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) |
---|
164 | zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) |
---|
165 | pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 |
---|
166 | pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 |
---|
167 | END DO |
---|
168 | END DO |
---|
169 | ELSE ! NO diffusion/viscosity |
---|
170 | CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' ) |
---|
171 | ENDIF |
---|
172 | ! ! deeper values (LAP and BLP cases) |
---|
173 | DO jk = 2, jpk |
---|
174 | pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk) |
---|
175 | pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk) |
---|
176 | END DO |
---|
177 | ! |
---|
178 | CASE( 'TRA' ) ! U- and V-points (approximation in zps-coord.) |
---|
179 | IF( cd_op == 'LAP' ) THEN ! laplacian operator |
---|
180 | IF(lwp) WRITE(numout,*) ' tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' |
---|
181 | za00 = pah0 / zd_max |
---|
182 | DO jj = 1, jpj |
---|
183 | DO ji = 1, jpi |
---|
184 | zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) |
---|
185 | zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) |
---|
186 | pah1(ji,jj,1) = za00 * zemax1 |
---|
187 | pah2(ji,jj,1) = za00 * zemax2 |
---|
188 | END DO |
---|
189 | END DO |
---|
190 | ELSEIF( cd_op == 'BLP' ) THEN ! bilaplacian operator (NB: square root of the coeff) |
---|
191 | IF(lwp) WRITE(numout,*) ' tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' |
---|
192 | za00 = pah0 / ( zd_max * zd_max * zd_max ) |
---|
193 | DO jj = 1, jpj |
---|
194 | DO ji = 1, jpi |
---|
195 | zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) |
---|
196 | zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) |
---|
197 | pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 |
---|
198 | pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 |
---|
199 | END DO |
---|
200 | END DO |
---|
201 | ELSE ! NO diffusion/viscosity |
---|
202 | CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' ) |
---|
203 | ENDIF |
---|
204 | ! ! deeper values (LAP and BLP cases) |
---|
205 | DO jk = 2, jpk |
---|
206 | pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk) |
---|
207 | pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk) |
---|
208 | END DO |
---|
209 | ! |
---|
210 | END SELECT |
---|
211 | ! |
---|
212 | END SUBROUTINE ldf_c2d |
---|
213 | |
---|
214 | !!====================================================================== |
---|
215 | END MODULE ldfc1d_c2d |
---|