1 | MODULE zdfmfc |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE zdfmfc *** |
---|
4 | !! Ocean physics: Mass-Flux scheme parameterization of Convection: |
---|
5 | !! Non-local transport for the convective ocean boundary |
---|
6 | !! layer. Subgrid-scale large eddies are represented by a |
---|
7 | !! mass-flux contribution (ln_zdfmfc = .TRUE.) |
---|
8 | !!====================================================================== |
---|
9 | !! History : NEMO ! |
---|
10 | !! 3.6 ! 2016-06 (H. Giordani, R. Bourdallé-Badie) Original code |
---|
11 | !! 4.2 ! 2020-12 (H. Giordani, R. Bourdallé-Badie) adapt to NEM04.2 |
---|
12 | !!---------------------------------------------------------------------- |
---|
13 | !!---------------------------------------------------------------------- |
---|
14 | !! tra_mfc : Compute the Mass Flux and trends of T/S |
---|
15 | !! diag_mfc : Modify diagonal of trazdf Matrix |
---|
16 | !! rhs_mfc : Modify RHS of trazdf Matrix |
---|
17 | !! zdf_mfc_init : initialization, namelist read, and parameters control |
---|
18 | !!---------------------------------------------------------------------- |
---|
19 | ! |
---|
20 | USE oce ! ocean dynamics and active tracers |
---|
21 | USE dom_oce ! ocean space and time domain |
---|
22 | USE domvvl ! ocean space and time domain : variable volume layer |
---|
23 | USE domzgr |
---|
24 | USE zdf_oce ! ocean vertical physics |
---|
25 | USE sbc_oce ! surface boundary condition: ocean |
---|
26 | USE phycst ! physical constants |
---|
27 | USE eosbn2 ! equation of state (eos routine) |
---|
28 | USE zdfmxl ! mixed layer |
---|
29 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
30 | USE lib_mpp ! MPP manager |
---|
31 | USE prtctl ! Print control |
---|
32 | USE in_out_manager ! I/O manager |
---|
33 | USE iom ! I/O manager library |
---|
34 | USE timing ! Timing |
---|
35 | USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) |
---|
36 | |
---|
37 | IMPLICIT NONE |
---|
38 | PRIVATE |
---|
39 | |
---|
40 | PUBLIC tra_mfc ! routine called in step module |
---|
41 | PUBLIC diag_mfc ! routine called in trazdf module |
---|
42 | PUBLIC rhs_mfc ! routine called in trazdf module |
---|
43 | PUBLIC zdf_mfc_init ! routine called in nemo module |
---|
44 | ! |
---|
45 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: edmfa, edmfb, edmfc !: diagonal term of the matrix. |
---|
46 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: edmftra !: y term for matrix inversion |
---|
47 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: edmfm !: y term for matrix inversion |
---|
48 | ! |
---|
49 | !! ** Namelist namzdf_edmf ** |
---|
50 | REAL(wp) :: rn_cemf ! entrain of T/S |
---|
51 | REAL(wp) :: rn_cwmf ! detrain of T/S |
---|
52 | REAL(wp) :: rn_cent ! entrain of the convective mass flux |
---|
53 | REAL(wp) :: rn_cdet ! detrain of the convective mass flux |
---|
54 | REAL(wp) :: rn_cap ! Factor of computation for convective area (negative => area constant) |
---|
55 | REAL(wp) :: App_max ! Maximum of the convective area |
---|
56 | LOGICAL, PUBLIC, SAVE :: ln_edmfuv !: EDMF flag for velocity ! |
---|
57 | ! |
---|
58 | !! * Substitutions |
---|
59 | # include "do_loop_substitute.h90" |
---|
60 | # include "domzgr_substitute.h90" |
---|
61 | # include "single_precision_substitute.h90" |
---|
62 | !!---------------------------------------------------------------------- |
---|
63 | !! NEMO/OCE 4.2 , NEMO Consortium (2018) |
---|
64 | !! $Id: zdfmfc.F90 13783 2020-20-02 15:30:22Z rbourdal $ |
---|
65 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
66 | !!---------------------------------------------------------------------- |
---|
67 | CONTAINS |
---|
68 | |
---|
69 | INTEGER FUNCTION zdf_mfc_alloc() |
---|
70 | !!---------------------------------------------------------------------- |
---|
71 | !! *** FUNCTION zdf_edmf_alloc *** |
---|
72 | !!---------------------------------------------------------------------- |
---|
73 | ALLOCATE( edmfa(jpi,jpj,jpk) , edmfb(jpi,jpj,jpk) , edmfc(jpi,jpj,jpk) & |
---|
74 | & , edmftra(jpi,jpj,jpk,2), edmfm(jpi,jpj,jpk) , STAT= zdf_mfc_alloc ) |
---|
75 | ! |
---|
76 | IF( lk_mpp ) CALL mpp_sum ( 'zdfmfc', zdf_mfc_alloc ) |
---|
77 | IF( zdf_mfc_alloc /= 0 ) CALL ctl_warn('zdf_mfc_alloc: failed to allocate arrays') |
---|
78 | END FUNCTION zdf_mfc_alloc |
---|
79 | |
---|
80 | |
---|
81 | SUBROUTINE tra_mfc( kt, Kmm, pts, Krhs ) |
---|
82 | !!---------------------------------------------------------------------- |
---|
83 | !! *** ROUTINE zdf_mfc *** |
---|
84 | !! |
---|
85 | !! ** Purpose : Compute a mass flux, depending on surface flux, over |
---|
86 | !! the instable part of the water column. |
---|
87 | !! |
---|
88 | !! ** Method : Compute surface instability and mix tracers until stable level |
---|
89 | !! |
---|
90 | !! |
---|
91 | !! ** Action : Compute convection plume and (ta,sa)-trends for trazdf (EDMF scheme) |
---|
92 | !! |
---|
93 | !! References : |
---|
94 | !! Giordani, Bourdallé-Badie and Madec JAMES 2020 |
---|
95 | !!---------------------------------------------------------------------- |
---|
96 | !!---------------------------------------------------------------------- |
---|
97 | INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices |
---|
98 | REAL(dp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation |
---|
99 | REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: ztsp ! T/S of the plume |
---|
100 | REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: ztse ! T/S at W point |
---|
101 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp ! |
---|
102 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp2 ! |
---|
103 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zapp ! |
---|
104 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zedmf ! |
---|
105 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zepsT, zepsW ! |
---|
106 | ! |
---|
107 | REAL(wp), DIMENSION(jpi,jpj) :: zustar, zustar2 ! |
---|
108 | REAL(wp), DIMENSION(jpi,jpj) :: zuws, zvws, zsws, zfnet ! |
---|
109 | REAL(wp), DIMENSION(jpi,jpj) :: zfbuo, zrautbm1 |
---|
110 | REAL(dp), DIMENSION(jpi,jpj) :: zrautb, zraupl |
---|
111 | REAL(wp), DIMENSION(jpi,jpj) :: zwpsurf ! |
---|
112 | REAL(wp), DIMENSION(jpi,jpj) :: zop0 , zsp0 ! |
---|
113 | REAL(wp), DIMENSION(jpi,jpj) :: zrwp_0, zrwp2_0 ! |
---|
114 | REAL(wp), DIMENSION(jpi,jpj) :: zapp0 ! |
---|
115 | REAL(wp), DIMENSION(jpi,jpj) :: zphp, zph, zphpm1, zphm1, zNHydro |
---|
116 | REAL(wp), DIMENSION(jpi,jpj) :: zhcmo ! |
---|
117 | ! |
---|
118 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zn2 ! N^2 |
---|
119 | REAL(wp), DIMENSION(jpi,jpj,2 ) :: zab, zabm1, zabp ! alpha and beta |
---|
120 | |
---|
121 | REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value |
---|
122 | |
---|
123 | REAL(wp) :: zrho, zrhop |
---|
124 | REAL(wp) :: zcnh, znum, zden, zcoef1, zcoef2 |
---|
125 | REAL(wp) :: zca, zcb, zcd, zrw, zxl, zcdet, zctre |
---|
126 | REAL(wp) :: zaw, zbw, zxw |
---|
127 | REAL(wp) :: alpha |
---|
128 | ! |
---|
129 | INTEGER, INTENT(in ) :: kt ! ocean time-step index ! |
---|
130 | ! |
---|
131 | INTEGER :: ji, jj, jk ! dummy loop arguments |
---|
132 | ! |
---|
133 | !------------------------------------------------------------------ |
---|
134 | ! Initialisation of coefficients |
---|
135 | !------------------------------------------------------------------ |
---|
136 | zca = 1._wp |
---|
137 | zcb = 1._wp |
---|
138 | zcd = 1._wp |
---|
139 | |
---|
140 | !------------------------------------------------------------------ |
---|
141 | ! Surface boundary condition |
---|
142 | !------------------------------------------------------------------ |
---|
143 | ! surface Stress |
---|
144 | !-------------------- |
---|
145 | zuws(:,:) = utau(:,:) * r1_rho0 |
---|
146 | zvws(:,:) = vtau(:,:) * r1_rho0 |
---|
147 | zustar2(:,:) = SQRT(zuws(:,:)*zuws(:,:)+zvws(:,:)*zvws(:,:)) |
---|
148 | zustar(:,:) = SQRT(zustar2(:,:)) |
---|
149 | |
---|
150 | ! Heat Flux |
---|
151 | !-------------------- |
---|
152 | zfnet(:,:) = qns(:,:) + qsr(:,:) |
---|
153 | zfnet(:,:) = zfnet(:,:) / (rho0 * rcp) |
---|
154 | |
---|
155 | ! Water Flux |
---|
156 | !--------------------- |
---|
157 | zsws(:,:) = emp(:,:) |
---|
158 | |
---|
159 | !------------------------------------------- |
---|
160 | ! Initialisation of prognostic variables |
---|
161 | !------------------------------------------- |
---|
162 | zrwp (:,:,:) = 0._wp ; zrwp2(:,:,:) = 0._wp ; zedmf(:,:,:) = 0._wp |
---|
163 | zph (:,:) = 0._wp ; zphm1(:,:) = 0._wp ; zphpm1(:,:) = 0._wp |
---|
164 | ztsp(:,:,:,:)= 0._wp |
---|
165 | |
---|
166 | ! Tracers inside plume (ztsp) and environment (ztse) |
---|
167 | ztsp(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1) |
---|
168 | ztsp(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) |
---|
169 | ztse(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1) |
---|
170 | ztse(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) |
---|
171 | |
---|
172 | CALL eos( ztse(:,:,1,:) , zrautb(:,:) ) |
---|
173 | CALL eos( ztsp(:,:,1,:) , zraupl(:,:) ) |
---|
174 | |
---|
175 | !------------------------------------------- |
---|
176 | ! Boundary Condition of Mass Flux (plume velo.; convective area, entrain/detrain) |
---|
177 | !------------------------------------------- |
---|
178 | zhcmo(:,:) = e3t(:,:,1,Kmm) |
---|
179 | zfbuo(:,:) = 0._wp |
---|
180 | WHERE ( ABS(zrautb(:,:)) > 1.e-20 ) zfbuo(:,:) = & |
---|
181 | & grav * ( 2.e-4_wp *zfnet(:,:) - 7.6E-4_wp*pts(:,:,1,jp_sal,Kmm)*zsws(:,:)/zrautb(:,:)) * zhcmo(:,:) |
---|
182 | |
---|
183 | zedmf(:,:,1) = -0.065_wp*(ABS(zfbuo(:,:)))**(1._wp/3._wp)*SIGN(1._wp,zfbuo(:,:)) |
---|
184 | zedmf(:,:,1) = MAX(0., zedmf(:,:,1)) |
---|
185 | |
---|
186 | zwpsurf(:,:) = 2._wp/3._wp*zustar(:,:) + 2._wp/3._wp*ABS(zfbuo(:,:))**(1._wp/3._wp) |
---|
187 | zwpsurf(:,:) = MAX(1.e-5_wp,zwpsurf(:,:)) |
---|
188 | zwpsurf(:,:) = MIN(1.,zwpsurf(:,:)) |
---|
189 | |
---|
190 | zapp(:,:,:) = App_max |
---|
191 | WHERE(zwpsurf .NE. 0.) zapp(:,:,1) = MIN(MAX(0.,zedmf(:,:,1)/zwpsurf(:,:)), App_max) |
---|
192 | |
---|
193 | zedmf(:,:,1) = 0._wp |
---|
194 | zrwp (:,:,1) = 0._wp |
---|
195 | zrwp2(:,:,1) = 0._wp |
---|
196 | zepsT(:,:,:) = 0.001_wp |
---|
197 | zepsW(:,:,:) = 0.001_wp |
---|
198 | |
---|
199 | |
---|
200 | !-------------------------------------------------------------- |
---|
201 | ! Compute plume properties |
---|
202 | ! In the same loop on vert. levels computation of: |
---|
203 | ! - Vertical velocity: zWp |
---|
204 | ! - Convective Area: zAp |
---|
205 | ! - Tracers properties inside the plume (if necessary): ztp |
---|
206 | !--------------------------------------------------------------- |
---|
207 | |
---|
208 | DO jk= 2, jpk |
---|
209 | |
---|
210 | ! Compute the buoyancy acceleration on T-points at jk-1 |
---|
211 | zrautbm1(:,:) = zrautb(:,:) |
---|
212 | CALL eos( CASTWP(pts (:,:,jk ,:,Kmm)) , zrautb(:,:) ) |
---|
213 | CALL eos( CASTWP(ztsp(:,:,jk-1,: )) , zraupl(:,:) ) |
---|
214 | |
---|
215 | zphm1(:,:) = zphm1(:,:) + grav * zrautbm1(:,:) * e3t(:,:,jk-1, Kmm) |
---|
216 | zphpm1(:,:) = zphpm1(:,:) + grav * zraupl(:,:) * e3t(:,:,jk-1, Kmm) |
---|
217 | zph(:,:) = zphm1(:,:) + grav * zrautb(:,:) * e3t(:,:,jk , Kmm) |
---|
218 | zph(:,:) = MAX( zph(:,:), zepsilon) |
---|
219 | |
---|
220 | WHERE(zrautbm1 .NE. 0.) zfbuo(:,:) = grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:) |
---|
221 | |
---|
222 | DO_2D( 0, 0, 0, 0 ) |
---|
223 | |
---|
224 | ! Compute Environment of Plume. Interpolation T/S (before time step) on W-points |
---|
225 | zrw = (gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm)) & |
---|
226 | & / (gdept(ji,jj,jk,Kmm) - gdept(ji,jj,jk-1,Kmm)) |
---|
227 | ztse(ji,jj,jk,:) = (pts(ji,jj,jk,:,Kmm) * zrw + pts(ji,jj,jk-1,:,Kmm)*(1._wp - zrw) )*tmask(ji,jj,jk) |
---|
228 | |
---|
229 | !--------------------------------------------------------------- |
---|
230 | ! Compute the vertical velocity on W-points |
---|
231 | !--------------------------------------------------------------- |
---|
232 | |
---|
233 | ! Non-hydrostatic pressure terms in the wp2 equation |
---|
234 | zcnh = 0.2_wp |
---|
235 | znum = 0.5_wp + zcnh - & |
---|
236 | (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) & |
---|
237 | *e3t(ji,jj,jk-1,Kmm)*0.5_wp |
---|
238 | zden = 0.5_wp + zcnh + & |
---|
239 | (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) & |
---|
240 | *e3t(ji,jj,jk-1,Kmm)*0.5_wp |
---|
241 | |
---|
242 | zcoef1 = zca*e3t(ji,jj,jk-1,Kmm) / zden |
---|
243 | zcoef2 = znum/zden |
---|
244 | |
---|
245 | ! compute wp2 |
---|
246 | zrwp2(ji,jj,jk) = zcoef1*zfbuo(ji,jj) & |
---|
247 | + zcoef2*zrwp2(ji,jj,jk-1) |
---|
248 | zrwp2(ji,jj,jk) = MAX ( zrwp2(ji,jj,jk)*wmask(ji,jj,jk) , 0.) |
---|
249 | zrwp (ji,jj,jk) = SQRT( zrwp2(ji,jj,jk) ) |
---|
250 | |
---|
251 | !---------------------------------------------------------------------------------- |
---|
252 | ! Compute convective area on W-point |
---|
253 | ! Compute vertical profil of the convective area with mass conservation hypothesis |
---|
254 | ! If rn_cap negative => constant value on the water column. |
---|
255 | !---------------------------------------------------------------------------------- |
---|
256 | IF( rn_cap .GT. 0. ) THEN |
---|
257 | |
---|
258 | zxw = MAX(zrwp(ji,jj,jk-1), zrwp(ji,jj,jk) ) |
---|
259 | IF( zxw > 0. ) THEN |
---|
260 | |
---|
261 | zxl = (zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk))/(e3t(ji,jj,jk-1,Kmm)*zxw) |
---|
262 | IF (zxl .LT. 0._wp) THEN |
---|
263 | zctre = -1.*rn_cap*zxl |
---|
264 | zcdet = 0._wp |
---|
265 | ELSE |
---|
266 | zctre = 0._wp |
---|
267 | zcdet = rn_cap*zxl |
---|
268 | END IF |
---|
269 | zapp(ji,jj,jk) = zapp(ji,jj,jk-1)* & |
---|
270 | & (1._wp + (zxl + zctre - zcdet )*e3t(ji,jj,jk-1,Kmm)) |
---|
271 | ELSE |
---|
272 | zapp(ji,jj,jk) = App_max |
---|
273 | END IF |
---|
274 | zapp(ji,jj,jk) = MIN( MAX(zapp(ji,jj,jk),0.), App_max) |
---|
275 | ELSE |
---|
276 | zapp(ji,jj,jk) = -1. * rn_cap |
---|
277 | END IF |
---|
278 | |
---|
279 | ! Compute Mass Flux on W-point |
---|
280 | zedmf(ji,jj,jk) = -zapp(ji,jj,jk) * zrwp(ji,jj,jk)* wmask(ji,jj,jk) |
---|
281 | |
---|
282 | ! Compute Entrainment coefficient |
---|
283 | IF(rn_cemf .GT. 0.) THEN |
---|
284 | zxw = 0.5_wp*(zrwp(ji,jj,jk-1)+ zrwp(ji,jj,jk) ) |
---|
285 | zepsT(ji,jj,jk) = 0.01_wp |
---|
286 | IF( zxw > 0. ) THEN |
---|
287 | zepsT(ji,jj,jk) = zepsT(ji,jj,jk) + & |
---|
288 | & ABS( zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk) ) & |
---|
289 | & / ( e3t(ji,jj,jk-1,Kmm) * zxw ) |
---|
290 | zepsT(ji,jj,jk) = zepsT(ji,jj,jk) * rn_cemf * wmask(ji,jj,jk) |
---|
291 | ENDIF |
---|
292 | ELSE |
---|
293 | zepsT(ji,jj,jk) = -rn_cemf |
---|
294 | ENDIF |
---|
295 | |
---|
296 | ! Compute the detrend coef for velocity (on W-point and not T-points, bug ???) |
---|
297 | IF(rn_cwmf .GT. 0.) THEN |
---|
298 | zepsW(ji,jj,jk) = rn_cwmf * zepsT(ji,jj,jk) |
---|
299 | ELSE |
---|
300 | zepsW(ji,jj,jk) = -rn_cwmf |
---|
301 | ENDIF |
---|
302 | |
---|
303 | !--------------------------------------------------------------- |
---|
304 | ! Compute the plume properties on T-points |
---|
305 | !--------------------------------------------------------------- |
---|
306 | IF(zrwp (ji,jj,jk) .LT. 1.e-12_wp .AND. zrwp (ji,jj,jk-1) .LT. 1.e-12_wp) THEN |
---|
307 | ztsp(ji,jj,jk-1,jp_tem) = pts(ji,jj,jk-1,jp_tem,Kmm) |
---|
308 | ztsp(ji,jj,jk-1,jp_sal) = pts(ji,jj,jk-1,jp_sal,Kmm) |
---|
309 | ENDIF |
---|
310 | |
---|
311 | zcoef1 = (1._wp-zepsT(ji,jj,jk)*(1._wp-zrw)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk ) ) & |
---|
312 | & / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk) ) |
---|
313 | ! |
---|
314 | zcoef2 = zepsT(ji,jj,jk)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk) & |
---|
315 | & / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk)) |
---|
316 | ! |
---|
317 | ztsp(ji,jj,jk,jp_tem) = (zcoef1 * ztsp(ji,jj,jk-1,jp_tem) + & |
---|
318 | & zcoef2 * ztse(ji,jj,jk ,jp_tem) )*tmask(ji,jj,jk) |
---|
319 | ztsp(ji,jj,jk,jp_sal) = (zcoef1 * ztsp(ji,jj,jk-1,jp_sal) + & |
---|
320 | & zcoef2 * ztse(ji,jj,jk ,jp_sal) )*tmask(ji,jj,jk) |
---|
321 | |
---|
322 | END_2D |
---|
323 | END DO ! end of loop on jpk |
---|
324 | |
---|
325 | ! Compute Mass Flux on T-point |
---|
326 | DO jk=1,jpk-1 |
---|
327 | edmfm(:,:,jk) = (zedmf(:,:,jk+1) + zedmf(:,:,jk) )*0.5_wp |
---|
328 | END DO |
---|
329 | edmfm(:,:,jpk) = zedmf(:,:,jpk) |
---|
330 | |
---|
331 | ! Save variable (on T point) |
---|
332 | CALL iom_put( "mf_Tp" , ztsp(:,:,:,jp_tem) ) ! Save plume temperature |
---|
333 | CALL iom_put( "mf_Sp" , ztsp(:,:,:,jp_sal) ) ! Save plume salinity |
---|
334 | CALL iom_put( "mf_mf" , edmfm(:,:,:) ) ! Save Mass Flux |
---|
335 | ! Save variable (on W point) |
---|
336 | CALL iom_put( "mf_wp" , zrwp (:,:,:) ) ! Save convective velocity in the plume |
---|
337 | CALL iom_put( "mf_app", zapp (:,:,:) ) ! Save convective area |
---|
338 | |
---|
339 | !================================================================================= |
---|
340 | ! Computation of a tridiagonal matrix and right hand side terms of the linear system |
---|
341 | !================================================================================= |
---|
342 | edmfa(:,:,:) = 0._wp |
---|
343 | edmfb(:,:,:) = 0._wp |
---|
344 | edmfc(:,:,:) = 0._wp |
---|
345 | edmftra(:,:,:,:) = 0._wp |
---|
346 | |
---|
347 | !--------------------------------------------------------------- |
---|
348 | ! Diagonal terms |
---|
349 | !--------------------------------------------------------------- |
---|
350 | DO jk=1,jpk-1 |
---|
351 | edmfa(:,:,jk) = 0._wp |
---|
352 | edmfb(:,:,jk) = -edmfm(:,:,jk ) / e3w(:,:,jk+1,Kmm) |
---|
353 | edmfc(:,:,jk) = edmfm(:,:,jk+1) / e3w(:,:,jk+1,Kmm) |
---|
354 | END DO |
---|
355 | edmfa(:,:,jpk) = -edmfm(:,:,jpk-1) / e3w(:,:,jpk,Kmm) |
---|
356 | edmfb(:,:,jpk) = edmfm(:,:,jpk ) / e3w(:,:,jpk,Kmm) |
---|
357 | edmfc(:,:,jpk) = 0._wp |
---|
358 | |
---|
359 | !--------------------------------------------------------------- |
---|
360 | ! right hand side term for Temperature |
---|
361 | !--------------------------------------------------------------- |
---|
362 | DO jk=1,jpk-1 |
---|
363 | edmftra(:,:,jk,1) = - edmfm(:,:,jk ) * ztsp(:,:,jk ,jp_tem) / e3w(:,:,jk+1,Kmm) & |
---|
364 | & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_tem) / e3w(:,:,jk+1,Kmm) |
---|
365 | END DO |
---|
366 | edmftra(:,:,jpk,1) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_tem) / e3w(:,:,jpk,Kmm) & |
---|
367 | & + edmfm(:,:,jpk ) * ztsp(:,:,jpk ,jp_tem) / e3w(:,:,jpk,Kmm) |
---|
368 | |
---|
369 | !--------------------------------------------------------------- |
---|
370 | ! Right hand side term for Salinity |
---|
371 | !--------------------------------------------------------------- |
---|
372 | DO jk=1,jpk-1 |
---|
373 | edmftra(:,:,jk,2) = - edmfm(:,:,jk ) * ztsp(:,:,jk ,jp_sal) / e3w(:,:,jk+1,Kmm) & |
---|
374 | & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_sal) / e3w(:,:,jk+1,Kmm) |
---|
375 | END DO |
---|
376 | edmftra(:,:,jpk,2) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_sal) / e3w(:,:,jpk,Kmm) & |
---|
377 | & + edmfm(:,:,jpk ) * ztsp(:,:,jpk ,jp_sal) / e3w(:,:,jpk,Kmm) |
---|
378 | ! |
---|
379 | ! |
---|
380 | CALL lbc_lnk( 'zdfmfc', edmfm,'T',1._wp, edmfa,'T',1._wp, edmfb,'T',1._wp, edmfc,'T',1._wp, edmftra(:,:,:,jp_tem),'T',1._wp, edmftra(:,:,:,jp_sal),'T',1._wp) |
---|
381 | ! |
---|
382 | END SUBROUTINE tra_mfc |
---|
383 | |
---|
384 | |
---|
385 | SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa ) |
---|
386 | |
---|
387 | REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: zdiagi, zdiagd, zdiags ! inout: tridaig. terms |
---|
388 | REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step |
---|
389 | INTEGER , INTENT(in ) :: Kaa ! ocean time level indices |
---|
390 | |
---|
391 | INTEGER :: ji, jj, jk ! dummy loop arguments |
---|
392 | |
---|
393 | DO_3D( 0, 0, 0, 0, 1, jpkm1 ) |
---|
394 | zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk) |
---|
395 | zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk) |
---|
396 | zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk) |
---|
397 | END_3D |
---|
398 | |
---|
399 | END SUBROUTINE diag_mfc |
---|
400 | |
---|
401 | SUBROUTINE rhs_mfc( zrhs, jjn ) |
---|
402 | |
---|
403 | REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: zrhs ! inout: rhs trend |
---|
404 | INTEGER , INTENT(in ) :: jjn ! tracer indices |
---|
405 | |
---|
406 | INTEGER :: ji, jj, jk ! dummy loop arguments |
---|
407 | |
---|
408 | DO_3D( 0, 0, 0, 0, 1, jpkm1 ) |
---|
409 | zrhs(ji,jj,jk) = zrhs(ji,jj,jk) + edmftra(ji,jj,jk,jjn) |
---|
410 | END_3D |
---|
411 | |
---|
412 | END SUBROUTINE rhs_mfc |
---|
413 | |
---|
414 | |
---|
415 | |
---|
416 | SUBROUTINE zdf_mfc_init |
---|
417 | !!---------------------------------------------------------------------- |
---|
418 | !! *** ROUTINE zdf_mfc_init *** |
---|
419 | !! |
---|
420 | !! ** Purpose : Initialization of the vertical eddy diffivity and |
---|
421 | !! mass flux |
---|
422 | !! |
---|
423 | !! ** Method : Read the namzdf_mfc namelist and check the parameters |
---|
424 | !! called at the first timestep (nit000) |
---|
425 | !! |
---|
426 | !! ** input : Namlist namzdf_mfc |
---|
427 | !! |
---|
428 | !! ** Action : Increase by 1 the nstop flag is setting problem encounter |
---|
429 | !! |
---|
430 | !!---------------------------------------------------------------------- |
---|
431 | ! |
---|
432 | INTEGER :: jk ! dummy loop indices |
---|
433 | INTEGER :: ios ! Local integer output status for namelist read |
---|
434 | REAL(wp):: zcr ! local scalar |
---|
435 | !! |
---|
436 | NAMELIST/namzdf_mfc/ ln_edmfuv, rn_cemf, rn_cwmf, rn_cent, rn_cdet, rn_cap, App_max |
---|
437 | !!---------------------------------------------------------- |
---|
438 | ! |
---|
439 | ! |
---|
440 | ! REWIND( numnam_ref ) ! Namelist namzdf_mfc in reference namelist : Vertical eddy diffivity mass flux |
---|
441 | READ ( numnam_ref, namzdf_mfc, IOSTAT = ios, ERR = 901) |
---|
442 | 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in reference namelist' ) |
---|
443 | |
---|
444 | ! REWIND( numnam_cfg ) ! Namelist namzdf_mfc in configuration namelist : Vertical eddy diffivity mass flux |
---|
445 | READ ( numnam_cfg, namzdf_mfc, IOSTAT = ios, ERR = 902 ) |
---|
446 | 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in configuration namelist' ) |
---|
447 | IF(lwm) WRITE ( numond, namzdf_mfc ) |
---|
448 | |
---|
449 | IF(lwp) THEN !* Control print |
---|
450 | WRITE(numout,*) |
---|
451 | WRITE(numout,*) 'zdf_mfc_init' |
---|
452 | WRITE(numout,*) '~~~~~~~~~~~~~' |
---|
453 | WRITE(numout,*) ' Namelist namzdf_mfc : set eddy diffusivity Mass Flux Convection' |
---|
454 | WRITE(numout,*) ' Apply mass flux on velocities (Not yet avail.) ln_edmfuv = ', ln_edmfuv |
---|
455 | WRITE(numout,*) ' Coeff for entrain/detrain T/S of plume (Neg => cte) rn_cemf = ', rn_cemf |
---|
456 | WRITE(numout,*) ' Coeff for entrain/detrain Wp of plume (Neg => cte) rn_cwmf = ', rn_cwmf |
---|
457 | WRITE(numout,*) ' Coeff for entrain/detrain area of plume rn_cap = ', rn_cap |
---|
458 | WRITE(numout,*) ' Coeff for entrain area of plume rn_cent = ', rn_cent |
---|
459 | WRITE(numout,*) ' Coeff for detrain area of plume rn_cdet = ', rn_cdet |
---|
460 | WRITE(numout,*) ' Max convective area App_max = ', App_max |
---|
461 | ENDIF |
---|
462 | !* allocate edmf arrays |
---|
463 | IF( zdf_mfc_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_edmf_init : unable to allocate arrays' ) |
---|
464 | edmfa(:,:,:) = 0._wp |
---|
465 | edmfb(:,:,:) = 0._wp |
---|
466 | edmfc(:,:,:) = 0._wp |
---|
467 | edmftra(:,:,:,:) = 0._wp |
---|
468 | ! |
---|
469 | END SUBROUTINE zdf_mfc_init |
---|
470 | |
---|
471 | !!====================================================================== |
---|
472 | |
---|
473 | !!====================================================================== |
---|
474 | END MODULE zdfmfc |
---|
475 | |
---|
476 | |
---|
477 | |
---|