New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfmfc.F90 in NEMO/trunk/src/OCE/ZDF – NEMO

source: NEMO/trunk/src/OCE/ZDF/zdfmfc.F90 @ 14985

Last change on this file since 14985 was 14834, checked in by hadcv, 3 years ago

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