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/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/ZDF/zdfmfc.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago

Merge trunk -r14642:HEAD
File size: 22.2 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#  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   !!----------------------------------------------------------------------
67CONTAINS
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)
442901   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 )
446902   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   !!======================================================================
474END MODULE zdfmfc
475
476
477
Note: See TracBrowser for help on using the repository browser.