source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90 @ 11305

Last change on this file since 11305 was 11305, checked in by smasson, 21 months ago

dev_r11265_ABL : first implementation of ABL, see #2131

File size: 43.8 KB
Line 
1MODULE ablmod
2   !!======================================================================
3   !!                       ***  MODULE  ablmod  ***
4   !! Surface module :  ABL computation to provide atmospheric data
5   !!                   for surface fluxes computation
6   !!======================================================================
7   !! History :  3.6  ! 2019-03  (F. Lemarié & G. Samson)  Original code
8   !!----------------------------------------------------------------------
9   
10   !!----------------------------------------------------------------------
11   !!   abl_stp       : ABL single column model
12   !!   abl_zdf_tke   : atmospheric vertical closure
13   !!----------------------------------------------------------------------
14   USE abl            ! ABL dynamics and tracers
15   USE par_abl        ! ABL constants
16
17   USE phycst         ! physical constants
18   USE dom_oce, ONLY  : tmask 
19   USE sbc_oce, ONLY  : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1
20   USE sbcblk         ! use some physical constants for flux computation
21   !
22   USE prtctl         ! Print control                    (prt_ctl routine)
23   USE iom            ! IOM library
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE timing         ! Timing
27
28   IMPLICIT NONE
29
30   PUBLIC   abl_stp   ! called by sbcabl.F90
31   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ustar2
32   !! * Substitutions
33#  include "vectopt_loop_substitute.h90"
34
35CONTAINS
36
37
38!===================================================================================================
39   SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &                            ! in
40              &            pu_dta, pv_dta, pt_dta, pq_dta, &
41              &            pslp_dta, pgu_dta, pgv_dta, &
42              &            pcd_du, psen, pevp,     &                            ! in/out
43              &            pwndm, ptaui, ptauj, ptaum )
44!---------------------------------------------------------------------------------------------------
45
46      !!---------------------------------------------------------------------
47      !!                    ***  ROUTINE abl_stp ***
48      !!
49      !! ** Purpose :   Time-integration of the ABL model
50      !!
51      !! ** Method  :   Compute atmospheric variables : vertical turbulence
52      !!                             + Coriolis term + newtonian relaxation
53      !!               
54      !! ** Action  : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
55      !!              - Advance tracers to time n+1 (Euler backward scheme)
56      !!              - Compute Coriolis term with forward-backward scheme (possibly with geostrophic guide)
57      !!              - Advance u,v to time n+1 (Euler backward scheme)
58      !!              - Apply newtonian relaxation on the dynamics and the tracers
59      !!              - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
60      !!
61      !!----------------------------------------------------------------------
62      INTEGER  , INTENT(in   )                   ::   kt         ! time step index
63      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psst       ! sea-surface temperature [Celsius]
64      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu       ! sea-surface u (U-point)
65      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point)
66      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq       ! sea-surface humidity
67      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pu_dta     ! large-scale windi
68      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pv_dta     ! large-scale windj
69      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgu_dta    ! large-scale hpgi
70      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgv_dta    ! large-scale hpgj
71      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pt_dta     ! large-scale pot. temp.
72      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pq_dta     ! large-scale humidity
73      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp_dta   ! sea-level pressure
74      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du     ! Cd x Du (T-point)
75      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   psen       ! Ch x Du
76      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pevp       ! Ce x Du
77      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd||     
78      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui      ! taux
79      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj      ! tauy
80      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau||
81      !
82      REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zrhoa, zwnd_i, zwnd_j
83      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   zFC
84      REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF   
85      REAL(wp), DIMENSION(1:jpi,  jptq  )        ::   zBC
86      !
87      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a
88      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_b
89      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_c
90      !
91      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka )   ::   z_cft  !--FL--to be removed after the test phase
92      !
93      INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
94      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1
95      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv
96      LOGICAL             ::   ln_old_coriolis  = .FALSE.              ! possibility to switch off Coriolis term       
97      !
98      !!---------------------------------------------------------------------     
99      !
100      IF(lwp .AND. kt == nit000) THEN                  ! control print
101         WRITE(numout,*)
102         WRITE(numout,*) 'abl_stp : ABL time stepping'
103         WRITE(numout,*) '~~~~~~'
104      ENDIF 
105      !
106      IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
107      !! Compute ustar squared as Cd || Uatm-Uoce ||^2
108      !! needed for surface boundary condition of TKE
109      !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk)
110      DO jj = 1,jpj
111         DO ji = 1,jpi
112            ustar2(ji,jj) = pCd_du(ji,jj) * pwndm(ji,jj)
113         END DO
114      END DO 
115      !
116      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
117      !                            !  1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
118      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
119
120      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
121         
122      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
123      !                            !  2 *** Advance tracers to time n+1
124      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
125     
126      !-------------
127      DO jj = 1, jpj    ! outer loop             !--> tq_abl computed on (1:jpi) x (1:jpj)
128      !-------------   
129         ! Compute matrix elements for interior points
130         DO jk = 3, jpkam1
131            DO ji = 1, jpi   ! vector opt.
132               z_elem_a( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )          ! lower-diagonal
133               z_elem_c( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )          ! upper-diagonal       
134               z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )  !       diagonal           
135            END DO
136         END DO 
137         ! Boundary conditions 
138         DO ji = 1, jpi   ! vector opt.
139            ! Neumann at the bottom         
140            z_elem_a( ji,     2              ) = 0._wp
141            z_elem_c( ji,     2              ) = - rdt * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                         
142            ! Homogeneous Neumann at the top
143            z_elem_a( ji,     jpka           ) = - rdt * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
144            z_elem_c( ji,     jpka           ) = 0._wp
145            z_elem_b( ji,     jpka           ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
146         END DO           
147                   
148         DO jtra = 1,jptq  ! loop on active tracers
149               
150            DO jk = 3, jpkam1
151               DO ji = fs_2, fs_jpim1                 
152                  tq_abl  ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl  ( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side
153               END DO
154            END DO
155
156            IF(jtra == jp_ta) THEN
157               DO ji = 1,jpi  ! surface boundary condition for temperature             
158                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * psen(ji, jj) 
159                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   &
160                    &                               + rdt * psen(ji, jj) * ( psst(ji, jj) + rt0 )               
161                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
162               END DO
163            ELSE   
164               DO ji = 1,jpi  ! surface boundary condition for humidity             
165                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * pevp(ji, jj) 
166                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   &
167                    &                               + rdt * pevp(ji, jj) * pssq(ji, jj)               
168                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
169               END DO                                     
170            END IF
171            !!
172            !! Matrix inversion
173            !! ----------------------------------------------------------
174            DO ji = 1,jpi           
175               zcff                       =  1._wp / z_elem_b( ji, 2 )
176               zCF   (ji,   2           ) = - zcff * z_elem_c( ji, 2 )
177               tq_abl(ji,jj,2,nt_a,jtra) =    zcff * tq_abl(ji,jj,2,nt_a,jtra)
178            END DO             
179
180            DO jk = 3, jpka           
181               DO ji = 1,jpi           
182                  zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
183                  zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
184                  tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk  ,nt_a,jtra)   &
185                  &                - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) )
186               END DO
187            END DO
188!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...         
189            DO jk = jpkam1,2,-1           
190               DO ji = 1,jpi 
191                  tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) +    &
192                     &                        zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
193               END DO
194            END DO
195         
196         END DO   !<-- loop on tracers                   
197         !!
198      !-------------
199      END DO             ! end outer loop
200      !-------------   
201     
202     
203      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
204      !                            !  3 *** Compute Coriolis term with geostrophic guide
205      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
206      !-------------       
207      DO jk = 2, jpka    ! outer loop
208      !-------------
209         !             
210         ! Advance u_abl & v_abl to time n+1
211         DO jj = 1, jpj
212            DO ji = 1, jpi           
213               zcff = ( ff_t(ji,jj) * rdt )*( ff_t(ji,jj) * rdt )  ! (f dt)**2
214     
215               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  &
216                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    &
217                  &                 +  rdt * ff_t(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  &
218                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff)
219                 
220               v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  &
221                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   &
222                  &                 -  rdt * ff_t(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) &
223                  &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)               
224            END DO
225         END DO 
226         !                                   
227      !-------------
228      END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
229      !-------------               
230      !
231      IF( ln_geos_winds ) THEN
232         DO jj = 1, jpj    ! outer loop       
233            DO jk = 1, jpka
234               DO ji = 1, jpi 
235                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
236                     &                      - rdt * e3t_abl(jk) * ff_t(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
237                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   &
238                     &                      + rdt * e3t_abl(jk) * ff_t(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
239               END DO
240            END DO
241         END DO     
242      END IF 
243      !-------------           
244      !
245      IF( ln_hpgls_frc ) THEN
246         DO jj = 1, jpj    ! outer loop       
247            DO jk = 1, jpka
248               DO ji = 1, jpi 
249                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgu_dta(ji,jj,jk)
250                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgv_dta(ji,jj,jk) 
251               ENDDO
252            ENDDO
253         ENDDO 
254      END IF
255
256      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
257      !                            !  4 *** Advance u,v to time n+1
258      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
259      !
260      !  Vertical diffusion for u_abl     
261      !-------------
262      DO jj = 1, jpj    ! outer loop
263      !-------------   
264
265         DO jk = 3, jpkam1
266            DO ji = 1, jpi 
267               z_elem_a( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal
268               z_elem_c( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal               
269               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal
270            END DO
271         END DO 
272           
273         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)           
274            !++ Surface boundary condition
275            z_elem_a( ji,     2    ) = 0._wp
276            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
277            zcff  = rdt * pcd_du(ji, jj)
278            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff 
279            u_abl( ji, jj,    2, nt_a ) = u_abl( ji, jj,    2, nt_a ) + 0.5_wp * zcff * ( pssu(ji-1, jj) + pssu(ji,jj) ) 
280            !++ Top Neumann B.C.
281            !z_elem_a( ji,     jpka ) = - 0.5_wp * rdt * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )
282            !z_elem_c( ji,     jpka ) = 0._wp
283            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
284            !++ Top Dirichlet B.C.
285            z_elem_a( ji,     jpka )       = 0._wp
286            z_elem_c( ji,     jpka )       = 0._wp
287            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
288            u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)             
289         END DO 
290         !!
291         !! Matrix inversion
292         !! ----------------------------------------------------------
293         DO ji = 2, jpi         
294            zcff                 =   1._wp / z_elem_b( ji, 2 )
295            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 ) 
296            u_abl (ji,jj,2,nt_a) =    zcff * u_abl(ji,jj,2,nt_a)
297         END DO             
298
299         DO jk = 3, jpka           
300            DO ji = 2, jpi             
301               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
302               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
303               u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk  ,nt_a)   &
304               &          - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
305            END DO
306         END DO
307           
308         DO jk = jpkam1,2,-1           
309            DO ji = 2, jpi
310               u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a)
311            END DO
312         END DO
313                                         
314      !-------------
315      END DO             ! end outer loop
316      !-------------   
317
318      !
319      !  Vertical diffusion for v_abl     
320      !-------------
321      DO jj = 2, jpj   ! outer loop
322      !-------------   
323         !
324         DO jk = 3, jpkam1
325            DO ji = 1, jpi   
326               z_elem_a( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
327               z_elem_c( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal             
328               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal
329            END DO
330         END DO
331
332         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)         
333            !++ Surface boundary condition
334            z_elem_a( ji,     2    ) = 0._wp
335            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
336            zcff  = rdt * pcd_du(ji, jj)
337            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff 
338            v_abl( ji, jj,    2, nt_a ) =  v_abl( ji, jj, 2, nt_a ) + 0.5_wp * zcff * ( pssv(ji, jj) + pssv(ji, jj-1) ) 
339            !++ Top Neumann B.C.           
340            !z_elem_a( ji,     jpka ) = -rdt * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
341            !z_elem_c( ji,     jpka ) = 0._wp
342            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
343            !++ Top Dirichlet B.C.
344            z_elem_a( ji,     jpka )       = 0._wp
345            z_elem_c( ji,     jpka )       = 0._wp
346            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
347            v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)             
348         END DO 
349         !!
350         !! Matrix inversion
351         !! ----------------------------------------------------------
352         DO ji = 1, jpi             
353            zcff                 =  1._wp / z_elem_b( ji, 2 )
354            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       ) 
355            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a ) 
356         END DO             
357
358         DO jk = 3, jpka           
359            DO ji = 1, jpi                   
360               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
361               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
362               v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk  ,nt_a)   &
363               &          - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
364            END DO
365         END DO
366           
367         DO jk = jpkam1,2,-1           
368            DO ji = 1, jpi   
369               v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a)
370            END DO
371         END DO
372         !                                         
373      !-------------
374      END DO             ! end outer loop
375      !-------------   
376 
377      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
378      !                            !  5 *** Apply nudging on the dynamics and the tracers
379      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
380      z_cft(:,:,:) = 0._wp     
381     
382      IF( nn_dyn_restore > 0  ) THEN
383         !-------------
384         DO jk = 2, jpka    ! outer loop
385         !-------------       
386            DO jj = 2, jpj
387               DO ji = 2, jpi
388                  zcff1  = pblh( ji, jj )
389                  zsig   = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )                       
390                  zsig  =            MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
391                  zmsk  = msk_abl(ji,jj)
392                  zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   &
393                     &  + jp_alp1_dyn * zsig    + jp_alp0_dyn
394                  zcff  = (1._wp-zmsk) + zmsk * rdt * zcff2   ! zcff = 1 for masked points
395                 
396                  u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   &
397                     &                               + zcff   * pu_dta( ji, jj, jk       )                     
398                  v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   &
399                     &                               + zcff   * pv_dta( ji, jj, jk       )
400                END DO
401            END DO   
402         !-------------
403         END DO             ! end outer loop
404         !-------------               
405      END IF
406
407      !-------------
408      DO jk = 2, jpka    ! outer loop
409      !-------------       
410         DO jj = 1,jpj
411            DO ji = 1,jpi 
412               zcff1  = pblh( ji, jj )
413               zsig   = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
414               zsig   =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
415               zmsk   = msk_abl(ji,jj)
416               zcff2  = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   &
417                  &   + jp_alp1_tra * zsig    + jp_alp0_tra
418               zcff   = (1._wp-zmsk) + zmsk * rdt * zcff2   ! zcff = 1 for masked points     
419               tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   &
420                  &                                       + zcff   * pt_dta( ji, jj, jk )
421               
422               tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa )   &
423                  &                                       + zcff   * pq_dta( ji, jj, jk )
424               
425               z_cft( ji, jj, jk ) = zcff
426            END DO
427         END DO
428      !-------------
429      END DO             ! end outer loop
430      !-------------
431      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
432      !                            !  6 *** MPI exchanges 
433      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
434      !
435      CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), kfillmode = jpfillnothing  )
436      CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T', 1., kfillmode = jpfillnothing  )   ! ++++ this should not be needed...
437      !
438      CALL iom_put ( 'u_abl',  u_abl(:,:,2:jpka,nt_a      ))
439      CALL iom_put ( 'v_abl',  v_abl(:,:,2:jpka,nt_a      ))
440      CALL iom_put ( 't_abl', tq_abl(:,:,2:jpka,nt_a,jp_ta))
441      CALL iom_put ( 'q_abl', tq_abl(:,:,2:jpka,nt_a,jp_qa))
442!      CALL iom_put ( 'coeft',  z_cft(:,:,2:jpka           ))
443      CALL iom_put (  'pblh',   pblh(:,:                  ))
444      !
445      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
446      !                            !  7 *** Finalize flux computation
447      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
448
449      DO jj = 1, jpj
450         DO ji = 1, jpi
451            ztemp             = tq_abl  ( ji, jj, 2, nt_a, jp_ta ) 
452            zhumi             = tq_abl  ( ji, jj, 2, nt_a, jp_qa ) 
453            zcff              = pslp_dta( ji, jj ) /   &              !<-- At this point ztemp and zhumi should not be zero ...
454               &                        (  R_dry*ztemp * ( 1._wp + rctv0*zhumi )  )
455            psen ( ji, jj )   =      cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )
456            pevp ( ji, jj )   = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) )
457            zrhoa( ji, jj )   = zcff             
458         END DO
459      END DO
460     
461      DO jj = 2, jpj
462         DO ji = 2, jpi   ! vect. opt.
463            zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji  ,jj) + pssu(ji-1,jj) ) 
464            zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj  ) + pssv(ji,jj-1) ) 
465         END DO
466      END DO
467      !
468      CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
469      !
470      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
471      DO jj = 1, jpj
472         DO ji = 1, jpi     
473            zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   &
474               &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj)  ) * msk_abl(ji,jj)
475            zztmp         = zrhoa(ji,jj) * pcd_du(ji,jj)
476           
477            pwndm (ji,jj) =         zcff
478            ptaum (ji,jj) = zztmp * zcff
479            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
480            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
481         END DO
482      END DO   
483      ! ... utau, vtau at U- and V_points, resp.
484      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
485      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
486      DO jj = 2, jpjm1
487         DO ji = 2, jpim1
488            zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
489            zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
490            ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) )
491            zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
492            zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
493            ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) )
494         END DO
495      END DO
496      !
497      CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1., ptauj(:,:), 'V', -1. )
498
499      CALL iom_put( "taum_oce", ptaum )
500
501      IF(ln_ctl) THEN
502         CALL prt_ctl( tab2d_1=pwndm  , clinfo1=' abl_stp: wndm   : ' )
503         CALL prt_ctl( tab2d_1=ptaui  , clinfo1=' abl_stp: utau   : ' )
504         CALL prt_ctl( tab2d_2=ptauj  , clinfo2=          'vtau   : ' )
505      ENDIF
506      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
507      !                            !  8 *** Swap time indices for the next timestep
508      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
509      nt_n = 1 + MOD( kt  , 2)
510      nt_a = 1 + MOD( kt+1, 2)
511      !
512!---------------------------------------------------------------------------------------------------
513   END SUBROUTINE abl_stp
514!===================================================================================================
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532!===================================================================================================
533   SUBROUTINE abl_zdf_tke( )
534!---------------------------------------------------------------------------------------------------
535
536      !!----------------------------------------------------------------------
537      !!                   ***  ROUTINE abl_zdf_tke  ***
538      !!
539      !! ** Purpose :   Time-step Turbulente Kinetic Energy (TKE) equation
540      !!
541      !! ** Method  : - source term due to shear
542      !!              - source term due to stratification
543      !!              - resolution of the TKE equation by inverting
544      !!                a tridiagonal linear system
545      !!
546      !! ** Action  : - en : now turbulent kinetic energy)
547      !!              - avmu, avmv : production of TKE by shear at u and v-points
548      !!                (= Kz dz[Ub] * dz[Un] )
549      !! ---------------------------------------------------------------------
550      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn 
551      INTEGER, DIMENSION(1:jpi          )     ::   ikbl
552      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
553      REAL(wp)                                ::   zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2
554      REAL(wp)                                ::   zdU2,zdV2     
555      REAL(wp)                                ::   zwndi,zwndj
556      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2
557      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2
558      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF       
559      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a
560      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b
561      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_c
562      LOGICAL                                 ::   ln_Patankar    = .FALSE.
563      LOGICAL                                 ::   ln_dumpvar     = .FALSE.
564      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl 
565      !
566      tind  = nt_n
567      ziRic = 1._wp / rn_Ric
568      ! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
569      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
570      !                            !  Advance TKE equation to time n+1
571      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
572      !-------------
573      DO jj = 1, jpj    ! outer loop
574      !-------------
575         !
576         ! Compute vertical shear         
577         DO jk = 2, jpkam1
578            DO ji = 1,jpi 
579               zcff        = 1.0_wp / e3w_abl( jk )**2               
580               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2                           
581               zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
582               zsh2(ji,jk) = zdU+zdV
583            END DO
584         END DO   
585         !
586         ! Compute brunt-vaisala frequency
587         DO jk = 2, jpkam1
588            DO ji = 1,jpi 
589               zcff  = grav * itvref / e3w_abl( jk ) 
590               zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta)
591               zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        &
592                  &   - tq_abl( ji, jj, jk  , tind, jp_ta) * tq_abl( ji, jj, jk  , tind, jp_qa)
593               zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi)
594            END DO
595         END DO 
596         !
597         ! Terms for the tridiagonal problem
598         DO jk = 2, jpkam1
599            DO ji = 1,jpi 
600               zshear       =                 zsh2( ji,     jk )   ! zsh2 is already multiplied by Avm_abl at this point
601               zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation         
602               zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
603               
604               z_elem_a( ji,     jk )   = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk   ) ! lower-diagonal
605               z_elem_c( ji,     jk )   = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal         
606               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE
607                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
608                     &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal       
609                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * ( zbuoy + zshear ) )             ! right-hand-side
610               ELSE
611                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
612                     &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal   
613                     &                     - e3w_abl(jk) * rdt * zbuoy   
614                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt *  zshear )             ! right-hand-side                     
615               END IF
616            END DO
617         END DO 
618                           
619         DO ji = 1,jpi    ! vector opt.           
620            zesrf   =  MAX( 4.63_wp * ustar2(ji,jj), tke_min )           
621            zetop   = tke_min             
622            z_elem_a ( ji,     1    ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp                                                   
623            z_elem_a ( ji,     jpka ) = 0._wp; z_elem_c ( ji,     jpka ) = 0._wp; z_elem_b ( ji,     jpka ) = 1._wp           
624            tke_abl( ji, jj,    1, nt_a ) = zesrf
625            tke_abl( ji, jj, jpka, nt_a ) = zetop 
626            zbn2(ji,jj,   1) = zbn2( ji,jj,   2) 
627            zsh2(ji,      1) = zsh2( ji,      2)                                     
628            zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 
629            zsh2(ji,   jpka) = zsh2( ji  , jpkam1)
630         END DO       
631         !!
632         !! Matrix inversion
633         !! ----------------------------------------------------------
634         DO ji = 1,jpi
635            zcff                  =  1._wp / z_elem_b( ji, 1 )
636            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       ) 
637            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a ) 
638         END DO             
639
640         DO jk = 2, jpka           
641            DO ji = 1,jpi
642               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
643               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
644               tke_abl(ji,jj,jk,nt_a) =   zcff * ( tke_abl(ji,jj,jk  ,nt_a)   &
645               &          - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
646            END DO
647         END DO
648           
649         DO jk = jpkam1,1,-1           
650            DO ji = 1,jpi
651               tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a)
652            END DO
653         END DO
654         
655!!FL should not be needed because of Patankar procedure 
656         tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
657
658         !!
659         !! Diagnose PBL height
660         !! ----------------------------------------------------------
661         
662         
663         !                                                       
664         ! arrays zRH, zFC and zCF are available at this point
665         ! and zFC(:, 1 ) = 0.
666         ! diagnose PBL height based on zsh2 and zbn2
667         zFC (  :  ,1) = 0._wp
668         ikbl( 1:jpi ) = 0 
669         
670         DO jk = 2,jpka
671            DO ji = 1, jpi 
672               zcff  = ghw_abl( jk-1 )
673               zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
674               zcff  = ghw_abl( jk   )
675               zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*(                 &
676                               zcff2 * ( zsh2( ji, jk  ) - ziRic * zbn2( ji, jj, jk   ) &
677                           - rn_Cek  * ( ff_t( ji, jj  ) * ff_t( ji, jj ) ) ) &
678                             + zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
679                           - rn_Cek  * ( ff_t( ji, jj  ) * ff_t( ji, jj ) ) ) &
680                           &                                                 )
681               IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
682            END DO
683         END DO
684         !
685         ! finalize the computation of the PBL height
686         DO ji = 1, jpi
687            jk = ikbl(ji)
688            IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
689               pblh( ji, jj ) =  (  ghw_abl( jk-1 ) * zFC( ji, jk   )       &
690                  &              -  ghw_abl( jk   ) * zFC( ji, jk-1 )       &
691                  &              ) / ( zFC( ji, jk   ) - zFC( ji, jk-1 ) )
692            ELSE IF( jk==2 ) THEN
693               pblh( ji, jj ) = ghw_abl(2   )
694            ELSE
695               pblh( ji, jj ) = ghw_abl(jpka)
696            END IF
697         END DO
698         ! Optional : could add pblh smoothing if pblh is noisy horizontally ...
699         !
700      !-------------     
701      END DO     
702      !-------------
703      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
704      !                            !  Diagnostic mixing length computation
705      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
706      !
707      SELECT CASE ( nn_amxl )
708      !
709      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom         
710#   define zlup zRH     
711#   define zldw zFC 
712         DO jj = 1, jpj     ! outer loop
713            !
714            DO ji = 1, jpi
715               mxl_abl ( ji, jj,    1 )  = 0._wp
716               mxl_abl ( ji, jj, jpka )  = mxl_min
717               zldw( ji,        1 )  = 0._wp             
718               zlup( ji,     jpka )  = 0._wp
719            END DO
720            !
721            DO jk = 2, jpkam1   
722               DO ji = 1, jpi 
723                  zbuoy             = MAX( zbn2(ji, jj, jk), rsmall ) 
724                  mxl_abl( ji, jj, jk ) = MAX( mxl_min,  &
725                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   ) 
726               END DO
727            END DO   
728            !
729            ! Limit mxl
730            DO jk = jpkam1,1,-1   
731               DO ji = 1, jpi 
732                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) )
733               END DO
734            END DO   
735            !
736            DO jk = 2, jpka
737               DO ji = 1, jpi 
738                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )   
739               END DO
740            END DO 
741            !
742            DO jk = 1, jpka
743               DO ji = 1, jpi
744                  mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )   
745               END DO
746            END DO                         
747            !
748         END DO
749#   undef zlup     
750#   undef zldw 
751         !
752         !
753      CASE ( 1 )             ! length-scale computed as the distance to the PBL height
754         DO jj = 1,jpj      ! outer loop
755            !
756            DO ji = 1, jpi   ! vector opt.
757               zcff = 1._wp / pblh( ji, jj )     ! inverse of hbl             
758               DO jk = 1, jpka             
759                  zsig  = MIN( zcff * ghw_abl( jk ), 1. )   
760                  zcff1 = pblh( ji, jj )                 
761                  mxl_abl( ji, jj, jk ) =  mxl_min                           &
762                     &  +  zsig         * (  amx1*zcff1 + bmx1*mxl_min ) &
763                     &  +  zsig *  zsig * (  amx2*zcff1 + bmx2*mxl_min ) &
764                     &  +  zsig**3      * (  amx3*zcff1 + bmx3*mxl_min ) &
765                     &  +  zsig**4      * (  amx4*zcff1 + bmx4*mxl_min ) &
766                     &  +  zsig**5      * (  amx5*zcff1 + bmx5*mxl_min )
767               END DO
768            END DO 
769            !           
770         END DO           
771         !
772      CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale
773         !
774#   define zlup zRH     
775#   define zldw zFC           
776! zCF is used for matrix inversion
777!             
778       DO jj = 1, jpj      ! outer loop
779         
780         DO ji = 1, jpi           
781            zlup( ji,    1 ) = mxl_min
782            zldw( ji,    1 ) = mxl_min               
783            zlup( ji, jpka ) = mxl_min
784            zldw( ji, jpka ) = mxl_min               
785         END DO           
786         
787         DO jk = 2,jpka-1
788            DO ji = 1, jpi
789               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
790               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)       
791            END DO
792         END DO         
793         !!
794         !! BL89 search for lup
795         !! ----------------------------------------------------------           
796         DO jk=2,jpka-1 
797            !
798            DO ji = 1, jpi
799               zCF(ji,1:jpka) = 0._wp
800               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
801               ln_foundl(ji ) = .false.
802            END DO   
803            !           
804            DO jkup=jk+1,jpka-1
805               DO ji = 1, jpi
806                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
807                     &               ( zbn2(ji,jj,jkup  )*(ghw_abl(jkup  )-ghw_abl(jk)) &
808                     &               + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )
809                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
810                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
811                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)                   
812                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
813                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
814                     zlup(ji,jk) = zcff               
815                     ln_foundl(ji) = .true.
816                  END IF
817               END DO
818            END DO
819            !
820         END DO   
821         !!
822         !! BL89 search for ldwn
823         !! ----------------------------------------------------------         
824         DO jk=2,jpka-1         
825            !
826            DO ji = 1, jpi
827               zCF(ji,1:jpka) = 0._wp
828               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
829               ln_foundl(ji ) = .false.
830            END DO 
831            !   
832            DO jkdwn=jk-1,1,-1
833               DO ji = 1, jpi             
834                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
835                     &               * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
836                                       + zbn2(ji,jj,jkdwn  )*(ghw_abl(jk)-ghw_abl(jkdwn  )) )   
837                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
838                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
839                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )             
840                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
841                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
842                     zldw(ji,jk) = zcff 
843                     ln_foundl(ji) = .true.               
844                  END IF                                   
845               END DO
846            END DO
847            !     
848         END DO
849
850         DO jk = 1, jpka
851            DO ji = 1, jpi 
852               mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min ) 
853            END DO
854         END DO 
855
856      END DO
857#   undef zlup     
858#   undef zldw           
859         !
860      END SELECT     
861      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
862      !                            !  Finalize the computation of turbulent visc./diff.
863      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
864     
865      !-------------
866      DO jj = 1, jpj     ! outer loop
867      !-------------
868         DO jk = 1, jpka   
869            DO ji = 1, jpi  ! vector opt.
870               zcff              = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk )  &
871                  &                                       * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) ) 
872               zcff2             =  1. / ( 1. + zcff )   !<-- phi_z(z)
873               zcff              = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
874!!FL: MAX function probably useless because of the definition of mxl_min             
875               Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   )
876               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )                             
877            END DO
878         END DO
879      !-------------         
880      END DO     
881      !-------------
882
883!---------------------------------------------------------------------------------------------------
884   END SUBROUTINE abl_zdf_tke
885!===================================================================================================
886
887
888!!======================================================================
889END MODULE ablmod
Note: See TracBrowser for help on using the repository browser.