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 @ 14433

Last change on this file since 14433 was 14433, checked in by smasson, 3 years ago

trunk: merge dev_r14312_MPI_Interface into the trunk, #2598
File size: 22.1 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(jpi,jpj,jpk,2) ::   ztsp         ! T/S of the plume
99      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   ztse         ! T/S at W point
100      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp          !
101      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp2         !
102      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zapp          !
103      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zedmf         !
104      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zepsT, zepsW  !
105      !
106      REAL(wp), DIMENSION(jpi,jpj) :: zustar, zustar2   !
107      REAL(wp), DIMENSION(jpi,jpj) :: zuws, zvws, zsws, zfnet          !
108      REAL(wp), DIMENSION(jpi,jpj) :: zfbuo, zrautbm1, zrautb, zraupl
109      REAL(wp), DIMENSION(jpi,jpj) :: zwpsurf            !
110      REAL(wp), DIMENSION(jpi,jpj) :: zop0 , zsp0 !
111      REAL(wp), DIMENSION(jpi,jpj) :: zrwp_0, zrwp2_0  !
112      REAL(wp), DIMENSION(jpi,jpj) :: zapp0           !
113      REAL(wp), DIMENSION(jpi,jpj) :: zphp, zph, zphpm1, zphm1, zNHydro
114      REAL(wp), DIMENSION(jpi,jpj) :: zhcmo          !
115      !
116      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zn2    ! N^2
117      REAL(wp), DIMENSION(jpi,jpj,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      !------------------------------------------------------------------
139      ! Surface boundary condition
140      !------------------------------------------------------------------
141      ! surface Stress
142      !--------------------
143      zuws(:,:) = utau(:,:) * r1_rho0 
144      zvws(:,:) = vtau(:,:) * r1_rho0 
145      zustar2(:,:) = SQRT(zuws(:,:)*zuws(:,:)+zvws(:,:)*zvws(:,:))
146      zustar(:,:)  = SQRT(zustar2(:,:))
147
148      ! Heat Flux
149      !--------------------
150      zfnet(:,:) = qns(:,:) + qsr(:,:)
151      zfnet(:,:) = zfnet(:,:) / (rho0 * rcp)
152
153      ! Water Flux
154      !---------------------
155      zsws(:,:) = emp(:,:)
156
157      !-------------------------------------------
158      ! Initialisation of prognostic variables
159      !-------------------------------------------
160      zrwp (:,:,:) =  0._wp ; zrwp2(:,:,:) =  0._wp ; zedmf(:,:,:) =  0._wp
161      zph  (:,:)   =  0._wp ; zphm1(:,:)   =  0._wp ; zphpm1(:,:)  =  0._wp
162      ztsp(:,:,:,:)=  0._wp
163
164      ! Tracers inside plume (ztsp) and environment (ztse)
165      ztsp(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1)
166      ztsp(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
167      ztse(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1)
168      ztse(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
169
170      CALL eos( ztse(:,:,1,:) ,  zrautb(:,:) )
171      CALL eos( ztsp(:,:,1,:) ,  zraupl(:,:) )
172
173      !-------------------------------------------
174      ! Boundary Condition of Mass Flux (plume velo.; convective area, entrain/detrain)
175      !-------------------------------------------
176      zhcmo(:,:) = e3t(:,:,1,Kmm)
177      zfbuo(:,:)   = 0._wp
178      WHERE ( ABS(zrautb(:,:)) > 1.e-20 ) zfbuo(:,:)   =   &
179         &      grav * ( 2.e-4_wp *zfnet(:,:) - 7.6E-4_wp*pts(:,:,1,jp_sal,Kmm)*zsws(:,:)/zrautb(:,:)) * zhcmo(:,:)
180
181      zedmf(:,:,1) = -0.065_wp*(ABS(zfbuo(:,:)))**(1._wp/3._wp)*SIGN(1.,zfbuo(:,:))
182      zedmf(:,:,1) = MAX(0., zedmf(:,:,1))
183
184      zwpsurf(:,:) = 2._wp/3._wp*zustar(:,:) + 2._wp/3._wp*ABS(zfbuo(:,:))**(1._wp/3._wp)
185      zwpsurf(:,:) = MAX(1.e-5_wp,zwpsurf(:,:))
186      zwpsurf(:,:) = MIN(1.,zwpsurf(:,:))
187
188      zapp(:,:,:)  = App_max
189      WHERE(zwpsurf .NE. 0.) zapp(:,:,1)   = MIN(MAX(0.,zedmf(:,:,1)/zwpsurf(:,:)), App_max)
190
191      zedmf(:,:,1) = 0._wp 
192      zrwp (:,:,1) = 0._wp 
193      zrwp2(:,:,1) = 0._wp
194      zepsT(:,:,:) = 0.001_wp
195      zepsW(:,:,:) = 0.001_wp
196
197
198      !--------------------------------------------------------------
199      ! Compute plume properties
200      ! In the same loop on vert. levels computation of:
201      !    - Vertical velocity: zWp
202      !    - Convective Area: zAp
203      !    - Tracers properties inside the plume (if necessary): ztp
204      !---------------------------------------------------------------
205
206      DO jk= 2, jpk
207
208         ! Compute the buoyancy acceleration on T-points at jk-1
209         zrautbm1(:,:) = zrautb(:,:)
210         CALL eos( pts (:,:,jk  ,:,Kmm) ,  zrautb(:,:)   )
211         CALL eos( ztsp(:,:,jk-1,:    ) ,  zraupl(:,:)   )
212
213         zphm1(:,:)  = zphm1(:,:)  + grav * zrautbm1(:,:) * e3t(:,:,jk-1, Kmm)
214         zphpm1(:,:) = zphpm1(:,:) + grav * zraupl(:,:)   * e3t(:,:,jk-1, Kmm)
215         zph(:,:)    = zphm1(:,:)  + grav * zrautb(:,:)   * e3t(:,:,jk  , Kmm)
216         zph(:,:)    = MAX( zph(:,:), zepsilon)
217
218         WHERE(zrautbm1 .NE. 0.) zfbuo(:,:)  =  grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:)
219
220         DO_2D( 0, 0, 0, 0 )
221
222            ! Compute Environment of Plume. Interpolation T/S (before time step) on W-points
223            zrw              =  (gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm)) &
224               &              / (gdept(ji,jj,jk,Kmm) - gdept(ji,jj,jk-1,Kmm))
225            ztse(ji,jj,jk,:) = (pts(ji,jj,jk,:,Kmm) * zrw + pts(ji,jj,jk-1,:,Kmm)*(1._wp - zrw) )*tmask(ji,jj,jk)
226
227            !---------------------------------------------------------------
228            ! Compute the vertical velocity on W-points
229            !---------------------------------------------------------------
230
231            ! Non-hydrostatic pressure terms in the wp2 equation
232            zcnh = 0.2_wp 
233            znum = 0.5_wp  + zcnh - &
234                   (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) &
235                   *e3t(ji,jj,jk-1,Kmm)*0.5_wp   
236            zden = 0.5_wp + zcnh + &
237                   (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) &
238                   *e3t(ji,jj,jk-1,Kmm)*0.5_wp   
239
240            zcoef1 = zca*e3t(ji,jj,jk-1,Kmm) / zden
241            zcoef2 = znum/zden
242
243            ! compute wp2
244            zrwp2(ji,jj,jk) = zcoef1*zfbuo(ji,jj) &
245                            + zcoef2*zrwp2(ji,jj,jk-1)
246            zrwp2(ji,jj,jk) = MAX ( zrwp2(ji,jj,jk)*wmask(ji,jj,jk) , 0.)
247            zrwp (ji,jj,jk) = SQRT( zrwp2(ji,jj,jk) )
248
249            !----------------------------------------------------------------------------------
250            ! Compute convective area on W-point
251            ! Compute vertical profil of the convective area with mass conservation hypothesis
252            ! If rn_cap negative => constant value on the water column.
253            !----------------------------------------------------------------------------------
254            IF( rn_cap .GT. 0. ) THEN
255
256               zxw = MAX(zrwp(ji,jj,jk-1), zrwp(ji,jj,jk) )
257               IF( zxw > 0. ) THEN
258
259                  zxl = (zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk))/(e3t(ji,jj,jk-1,Kmm)*zxw)
260                  IF (zxl .LT. 0._wp) THEN
261                     zctre  = -1.*rn_cap*zxl 
262                     zcdet  =  0._wp
263                  ELSE
264                     zctre  =  0._wp
265                     zcdet  =  rn_cap*zxl 
266                  END IF
267                     zapp(ji,jj,jk) = zapp(ji,jj,jk-1)*     &
268                     &                (1._wp + (zxl + zctre - zcdet )*e3t(ji,jj,jk-1,Kmm))
269                  ELSE
270                     zapp(ji,jj,jk) = App_max
271                  END IF
272                  zapp(ji,jj,jk) = MIN( MAX(zapp(ji,jj,jk),0.), App_max)
273               ELSE
274                  zapp(ji,jj,jk) = -1. * rn_cap
275               END IF
276
277            ! Compute Mass Flux on W-point
278            zedmf(ji,jj,jk)   = -zapp(ji,jj,jk) * zrwp(ji,jj,jk)* wmask(ji,jj,jk)
279
280            ! Compute Entrainment coefficient
281            IF(rn_cemf .GT. 0.) THEN
282               zxw = 0.5_wp*(zrwp(ji,jj,jk-1)+ zrwp(ji,jj,jk) )
283               zepsT(ji,jj,jk)  =  0.01_wp
284               IF( zxw > 0.  ) THEN
285                  zepsT(ji,jj,jk)  =  zepsT(ji,jj,jk) +                       &
286                                   &  ABS( zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk) )  &
287                                   &  / ( e3t(ji,jj,jk-1,Kmm) * zxw )
288                  zepsT(ji,jj,jk)  = zepsT(ji,jj,jk) * rn_cemf * wmask(ji,jj,jk)
289               ENDIF
290            ELSE
291               zepsT(ji,jj,jk)  = -rn_cemf
292            ENDIF
293
294            ! Compute the detrend coef for velocity (on W-point and not T-points, bug ???)
295            IF(rn_cwmf .GT. 0.) THEN
296               zepsW(ji,jj,jk)  =  rn_cwmf * zepsT(ji,jj,jk)
297            ELSE
298               zepsW(ji,jj,jk)  = -rn_cwmf
299            ENDIF
300
301            !---------------------------------------------------------------
302            ! Compute the plume properties on T-points
303            !---------------------------------------------------------------
304            IF(zrwp (ji,jj,jk) .LT. 1.e-12_wp .AND. zrwp (ji,jj,jk-1) .LT. 1.e-12_wp) THEN
305               ztsp(ji,jj,jk-1,jp_tem) = pts(ji,jj,jk-1,jp_tem,Kmm)
306               ztsp(ji,jj,jk-1,jp_sal) = pts(ji,jj,jk-1,jp_sal,Kmm)
307            ENDIF
308
309            zcoef1 =  (1._wp-zepsT(ji,jj,jk)*(1._wp-zrw)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk ) ) &
310            &       / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk) )
311            !
312            zcoef2 =  zepsT(ji,jj,jk)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk)                       &
313            &       / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk))
314            !
315            ztsp(ji,jj,jk,jp_tem) = (zcoef1 * ztsp(ji,jj,jk-1,jp_tem) +  &
316            &                        zcoef2 * ztse(ji,jj,jk  ,jp_tem) )*tmask(ji,jj,jk)
317            ztsp(ji,jj,jk,jp_sal) = (zcoef1 * ztsp(ji,jj,jk-1,jp_sal) +  &
318            &                        zcoef2 * ztse(ji,jj,jk  ,jp_sal) )*tmask(ji,jj,jk)
319
320         END_2D 
321      END DO ! end of loop on jpk
322
323      ! Compute Mass Flux on T-point
324      DO jk=1,jpk-1
325         edmfm(:,:,jk) = (zedmf(:,:,jk+1)  + zedmf(:,:,jk) )*0.5_wp
326      END DO
327      edmfm(:,:,jpk) = zedmf(:,:,jpk) 
328
329      ! Save variable (on T point)
330      CALL iom_put( "mf_Tp" , ztsp(:,:,:,jp_tem) )  ! Save plume temperature
331      CALL iom_put( "mf_Sp" , ztsp(:,:,:,jp_sal) )  ! Save plume salinity
332      CALL iom_put( "mf_mf" , edmfm(:,:,:)       )  ! Save Mass Flux
333      ! Save variable (on W point)
334      CALL iom_put( "mf_wp" , zrwp (:,:,:)       )  ! Save convective velocity in the plume
335      CALL iom_put( "mf_app", zapp (:,:,:)       )  ! Save convective area
336
337      !=================================================================================
338      !  Computation of a tridiagonal matrix and right hand side terms of the linear system
339      !=================================================================================
340      edmfa(:,:,:)     = 0._wp
341      edmfb(:,:,:)     = 0._wp
342      edmfc(:,:,:)     = 0._wp
343      edmftra(:,:,:,:) = 0._wp
344
345      !---------------------------------------------------------------
346      ! Diagonal terms
347      !---------------------------------------------------------------
348      DO jk=1,jpk-1
349         edmfa(:,:,jk) =  0._wp
350         edmfb(:,:,jk) = -edmfm(:,:,jk  ) / e3w(:,:,jk+1,Kmm)
351         edmfc(:,:,jk) =  edmfm(:,:,jk+1) / e3w(:,:,jk+1,Kmm)
352      END DO
353      edmfa(:,:,jpk)   = -edmfm(:,:,jpk-1) / e3w(:,:,jpk,Kmm)
354      edmfb(:,:,jpk)   =  edmfm(:,:,jpk  ) / e3w(:,:,jpk,Kmm)
355      edmfc(:,:,jpk)   =  0._wp
356
357      !---------------------------------------------------------------
358      ! right hand side term for Temperature
359      !---------------------------------------------------------------
360      DO jk=1,jpk-1
361        edmftra(:,:,jk,1) = - edmfm(:,:,jk  ) * ztsp(:,:,jk  ,jp_tem) / e3w(:,:,jk+1,Kmm) &
362                          & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_tem) / e3w(:,:,jk+1,Kmm)
363      END DO
364      edmftra(:,:,jpk,1) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_tem) / e3w(:,:,jpk,Kmm) &
365                         & + edmfm(:,:,jpk  ) * ztsp(:,:,jpk  ,jp_tem) / e3w(:,:,jpk,Kmm)
366                       
367      !---------------------------------------------------------------
368      ! Right hand side term for Salinity
369      !---------------------------------------------------------------
370      DO jk=1,jpk-1
371         edmftra(:,:,jk,2) =  - edmfm(:,:,jk  ) * ztsp(:,:,jk  ,jp_sal) / e3w(:,:,jk+1,Kmm) &
372                           &  + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_sal) / e3w(:,:,jk+1,Kmm)
373      END DO
374      edmftra(:,:,jpk,2) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_sal) / e3w(:,:,jpk,Kmm) &
375                         & + edmfm(:,:,jpk  ) * ztsp(:,:,jpk  ,jp_sal) / e3w(:,:,jpk,Kmm)
376      !
377      !
378      CALL lbc_lnk( 'zdfmfc', edmfm,'T',1., edmfa,'T',1., edmfb,'T',1., edmfc,'T',1., edmftra(:,:,:,jp_tem),'T',1., edmftra(:,:,:,jp_sal),'T',1.)
379      !
380   END SUBROUTINE tra_mfc
381
382   
383   SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa )
384
385      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  zdiagi, zdiagd, zdiags  ! inout: tridaig. terms
386      REAL(wp)                        , INTENT(in   ) ::   p2dt                   ! tracer time-step
387      INTEGER                         , INTENT(in   ) ::   Kaa                    ! ocean time level indices
388
389      INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
390
391         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
392            zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk)
393            zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk)
394            zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk)
395         END_3D
396
397   END SUBROUTINE diag_mfc
398
399   SUBROUTINE rhs_mfc( zrhs, jjn )
400
401      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   zrhs                   ! inout: rhs trend
402      INTEGER                         , INTENT(in   ) ::   jjn                    ! tracer indices
403
404      INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
405
406      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
407         zrhs(ji,jj,jk) = zrhs(ji,jj,jk) + edmftra(ji,jj,jk,jjn)
408      END_3D
409
410   END SUBROUTINE rhs_mfc
411
412
413
414   SUBROUTINE zdf_mfc_init
415      !!----------------------------------------------------------------------
416      !!                  ***  ROUTINE zdf_mfc_init  ***
417      !!                   
418      !! ** Purpose :   Initialization of the vertical eddy diffivity and
419      !!      mass flux
420      !!
421      !! ** Method  :   Read the namzdf_mfc namelist and check the parameters
422      !!      called at the first timestep (nit000)
423      !!
424      !! ** input   :   Namlist namzdf_mfc
425      !!
426      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
427      !!
428      !!----------------------------------------------------------------------
429      !
430      INTEGER ::   jk    ! dummy loop indices
431      INTEGER ::   ios   ! Local integer output status for namelist read
432      REAL(wp)::   zcr   ! local scalar
433      !!
434      NAMELIST/namzdf_mfc/ ln_edmfuv, rn_cemf, rn_cwmf, rn_cent, rn_cdet, rn_cap, App_max
435      !!----------------------------------------------------------
436      !
437      !
438!      REWIND( numnam_ref )              ! Namelist namzdf_mfc in reference namelist : Vertical eddy diffivity mass flux
439      READ  ( numnam_ref, namzdf_mfc, IOSTAT = ios, ERR = 901)
440901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in reference namelist' )
441
442!      REWIND( numnam_cfg )              ! Namelist namzdf_mfc in configuration namelist : Vertical eddy diffivity mass flux
443      READ  ( numnam_cfg, namzdf_mfc, IOSTAT = ios, ERR = 902 )
444902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in configuration namelist' )
445      IF(lwm) WRITE ( numond, namzdf_mfc )
446
447      IF(lwp) THEN                     !* Control print
448         WRITE(numout,*)
449         WRITE(numout,*) 'zdf_mfc_init'
450         WRITE(numout,*) '~~~~~~~~~~~~~'
451         WRITE(numout,*) '   Namelist namzdf_mfc : set eddy diffusivity Mass Flux Convection'
452         WRITE(numout,*) '   Apply mass flux on velocities (Not yet avail.)     ln_edmfuv = ', ln_edmfuv
453         WRITE(numout,*) '   Coeff for entrain/detrain T/S of plume (Neg => cte) rn_cemf  = ', rn_cemf
454         WRITE(numout,*) '   Coeff for entrain/detrain Wp of plume  (Neg => cte) rn_cwmf  = ', rn_cwmf
455         WRITE(numout,*) '   Coeff for entrain/detrain area of plume             rn_cap   = ', rn_cap
456         WRITE(numout,*) '   Coeff for entrain area of plume                     rn_cent  = ', rn_cent
457         WRITE(numout,*) '   Coeff for detrain area of plume                     rn_cdet  = ', rn_cdet
458         WRITE(numout,*) '   Max convective area                                 App_max  = ', App_max
459       ENDIF
460                                     !* allocate edmf arrays
461      IF( zdf_mfc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_edmf_init : unable to allocate arrays' )
462      edmfa(:,:,:)     = 0._wp
463      edmfb(:,:,:)     = 0._wp
464      edmfc(:,:,:)     = 0._wp
465      edmftra(:,:,:,:) = 0._wp
466      !
467   END SUBROUTINE zdf_mfc_init
468
469   !!======================================================================
470 
471   !!======================================================================
472END MODULE zdfmfc
473
474
475
Note: See TracBrowser for help on using the repository browser.