source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/ablmod.F90 @ 13217

Last change on this file since 13217 was 13217, checked in by smasson, 3 months ago

better e3: update with trunk@13215 see #2385

File size: 64.4 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, rhoa
20   USE sbcblk         ! use rn_efac, cdn_oce
21   USE sbcblk_phy     ! use some physical constants for flux computation
22   !
23   USE prtctl         ! Print control                    (prt_ctl routine)
24   USE iom            ! IOM library
25   USE in_out_manager ! I/O manager
26   USE lib_mpp        ! MPP library
27   USE timing         ! Timing
28
29   IMPLICIT NONE
30
31   PUBLIC   abl_stp   ! called by sbcabl.F90
32   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ustar2, zrough
33   !! * Substitutions
34#  include "do_loop_substitute.h90"
35
36CONTAINS
37
38
39!===================================================================================================
40   SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq,            &     ! in
41              &            pu_dta, pv_dta, pt_dta, pq_dta,    &
42              &            pslp_dta, pgu_dta, pgv_dta,        &
43              &            pcd_du, psen, pevp,                &     ! in/out
44              &            pwndm, ptaui, ptauj, ptaum         &
45#if defined key_si3
46              &          , ptm_su, pssu_ice, pssv_ice         &
47              &          , pssq_ice, pcd_du_ice, psen_ice     &
48              &          , pevp_ice, pwndm_ice, pfrac_oce     &
49              &          , ptaui_ice, ptauj_ice               &
50#endif
51              &      )
52!---------------------------------------------------------------------------------------------------
53
54      !!---------------------------------------------------------------------
55      !!                    ***  ROUTINE abl_stp ***
56      !!
57      !! ** Purpose :   Time-integration of the ABL model
58      !!
59      !! ** Method  :   Compute atmospheric variables : vertical turbulence
60      !!                             + Coriolis term + newtonian relaxation
61      !!
62      !! ** Action  : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
63      !!              - Advance tracers to time n+1 (Euler backward scheme)
64      !!              - Compute Coriolis term with forward-backward scheme (possibly with geostrophic guide)
65      !!              - Advance u,v to time n+1 (Euler backward scheme)
66      !!              - Apply newtonian relaxation on the dynamics and the tracers
67      !!              - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
68      !!
69      !!----------------------------------------------------------------------
70      INTEGER  , INTENT(in   )                   ::   kt         ! time step index
71      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psst       ! sea-surface temperature [Celsius]
72      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu       ! sea-surface u (U-point)
73      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point)
74      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq       ! sea-surface humidity
75      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pu_dta     ! large-scale windi
76      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pv_dta     ! large-scale windj
77      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgu_dta    ! large-scale hpgi
78      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgv_dta    ! large-scale hpgj
79      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pt_dta     ! large-scale pot. temp.
80      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pq_dta     ! large-scale humidity
81      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp_dta   ! sea-level pressure
82      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du     ! Cd x Du (T-point)
83      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   psen       ! Ch x Du
84      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pevp       ! Ce x Du
85      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd||
86      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui      ! taux
87      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj      ! tauy
88      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau||
89      !
90#if defined key_si3
91      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptm_su       ! ice-surface temperature [K]
92      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu_ice     ! ice-surface u (U-point)
93      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point)
94      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity
95      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du_ice   ! Cd x Du over ice (T-point)
96      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psen_ice     ! Ch x Du over ice (T-point)
97      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pevp_ice     ! Ce x Du over ice (T-point)
98      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice||
99      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pfrac_oce    ! ocean fraction
100      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui_ice    ! ice-surface taux stress (U-point)
101      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point)
102#endif
103      !
104      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zwnd_i, zwnd_j
105      REAL(wp), DIMENSION(1:jpi      ,2:jpka)    ::   zCF
106      !
107      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_a
108      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_b
109      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_c
110      !
111      INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
112      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
113      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
114      LOGICAL             ::   SemiImp_Cor = .TRUE.
115      !
116      !!---------------------------------------------------------------------
117      !
118      IF(lwp .AND. kt == nit000) THEN                  ! control print
119         WRITE(numout,*)
120         WRITE(numout,*) 'abl_stp : ABL time stepping'
121         WRITE(numout,*) '~~~~~~'
122      ENDIF
123      !
124      IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
125      IF( kt == nit000 ) ALLOCATE ( zrough( 1:jpi, 1:jpj ) )
126      !! Compute ustar squared as Cd || Uatm-Uoce ||^2
127      !! needed for surface boundary condition of TKE
128      !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk)
129      DO_2D_11_11
130         zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj)
131#if defined key_si3
132         zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj)
133         ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice
134#else
135         ustar2(ji,jj) = zzoce
136#endif
137         zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(ji,jj), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce
138      END_2D
139      !
140      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
141      !                            !  1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
142      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
143
144      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
145
146      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
147      !                            !  2 *** Advance tracers to time n+1
148      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
149
150      !-------------
151      DO jj = 1, jpj    ! outer loop             !--> tq_abl computed on (1:jpi) x (1:jpj)
152      !-------------
153         ! Compute matrix elements for interior points
154         DO jk = 3, jpkam1
155            DO ji = 1, jpi   ! vector opt.
156               z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
157               z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal
158               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal
159            END DO
160         END DO
161         ! Boundary conditions
162         DO ji = 1, jpi   ! vector opt.
163            ! Neumann at the bottom
164            z_elem_a( ji,    2 ) = 0._wp
165            z_elem_c( ji,    2 ) = - rDt_abl * Avt_abl( ji, jj,    2 ) / e3w_abl(    2 )
166            ! Homogeneous Neumann at the top
167            z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )
168            z_elem_c( ji, jpka ) = 0._wp
169            z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
170         END DO
171
172         DO jtra = 1,jptq  ! loop on active tracers
173
174            DO jk = 3, jpkam1
175               !DO ji = 2, jpim1
176               DO ji = 1,jpi  !!GS: to be checked if needed
177                  tq_abl( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side
178               END DO
179            END DO
180
181            IF(jtra == jp_ta) THEN
182               DO ji = 1,jpi  ! surface boundary condition for temperature
183                  zztmp1 = psen(ji, jj)
184                  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )
185#if defined key_si3
186                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)
187                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
188#endif
189                  z_elem_b( ji,        2             ) = e3t_abl(    2 ) - z_elem_c( ji,        2             ) + rDt_abl * zztmp1
190                  tq_abl  ( ji, jj,    2, nt_a, jtra ) = e3t_abl(    2 ) * tq_abl  ( ji, jj,    2, nt_n, jtra ) + rDt_abl * zztmp2
191                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
192               END DO
193            ELSE ! jp_qa
194               DO ji = 1,jpi  ! surface boundary condition for humidity
195                  zztmp1 = pevp(ji, jj)
196                  zztmp2 = pevp(ji, jj) * pssq(ji, jj)
197#if defined key_si3
198                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj)
199                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)
200#endif
201                  z_elem_b( ji,     2                ) = e3t_abl(    2 ) - z_elem_c( ji,        2             ) + rDt_abl * zztmp1
202                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl(    2 ) * tq_abl  ( ji, jj,    2, nt_n, jtra ) + rDt_abl * zztmp2
203                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
204               END DO
205            END IF
206            !!
207            !! Matrix inversion
208            !! ----------------------------------------------------------
209            DO ji = 1,jpi
210               zcff                            =  1._wp / z_elem_b( ji, 2 )
211               zCF   ( ji,     2             ) = - zcff * z_elem_c( ji, 2 )
212               tq_abl( ji, jj, 2, nt_a, jtra ) =   zcff * tq_abl( ji, jj, 2, nt_a, jtra )
213            END DO
214
215            DO jk = 3, jpka
216               DO ji = 1,jpi
217                  zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF( ji, jk-1 ) )
218                  zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
219                  tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk  ,nt_a,jtra)   &
220                     &           - z_elem_a(ji, jk) *   tq_abl(ji,jj,jk-1,nt_a,jtra) )
221               END DO
222            END DO
223            !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...
224            DO jk = jpkam1,2,-1
225               DO ji = 1,jpi
226                  tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) +    &
227                     &                        zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
228               END DO
229            END DO
230
231         END DO   !<-- loop on tracers
232         !!
233      !-------------
234      END DO             ! end outer loop
235      !-------------
236
237      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
238      !                            !  3 *** Compute Coriolis term with geostrophic guide
239      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
240      IF( SemiImp_Cor ) THEN
241
242         !-------------
243         DO jk = 2, jpka    ! outer loop
244         !-------------
245            !
246            ! Advance u_abl & v_abl to time n+1
247            DO_2D_11_11
248               zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2
249
250               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(                                          &
251                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * u_abl( ji, jj, jk, nt_n )    &
252                  &                     + rDt_abl * fft_abl(ji, jj) * v_abl( ji, jj, jk, nt_n ) )  &
253                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff)
254   
255               v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(                                         &
256                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * v_abl( ji, jj, jk, nt_n )    &
257                  &                     - rDt_abl * fft_abl(ji, jj) * u_abl( ji, jj, jk, nt_n ) )  &
258                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff)
259            END_2D
260            !
261         !-------------
262         END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
263         !-------------
264         !
265         IF( ln_geos_winds ) THEN
266            DO jj = 1, jpj    ! outer loop
267               DO jk = 1, jpka
268                  DO ji = 1, jpi
269                     u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
270                        &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
271                     v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   &
272                        &                      + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
273                  END DO
274               END DO
275            END DO
276         END IF
277         !
278         IF( ln_hpgls_frc ) THEN
279            DO jj = 1, jpj    ! outer loop
280               DO jk = 1, jpka
281                  DO ji = 1, jpi
282                     u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk)
283                     v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)
284                  ENDDO
285               ENDDO
286            ENDDO
287         END IF
288
289      ELSE ! SemiImp_Cor = .FALSE.
290
291         IF( ln_geos_winds ) THEN
292
293            !-------------
294            DO jk = 2, jpka    ! outer loop
295            !-------------
296               !
297               IF( MOD( kt, 2 ) == 0 ) then
298                  ! Advance u_abl & v_abl to time n+1
299                  DO jj = 1, jpj
300                     DO ji = 1, jpi
301                        zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_n ) - pgv_dta(ji  ,jj  ,jk)  )
302                        u_abl( ji, jj, jk, nt_a ) =                u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff
303                        zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_a ) - pgu_dta(ji  ,jj  ,jk)  )
304                        v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff )
305                        u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  u_abl( ji, jj, jk, nt_a )
306                     END DO
307                  END DO
308               ELSE
309                  DO jj = 1, jpj
310                     DO ji = 1, jpi
311                        zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_n ) - pgu_dta(ji  ,jj  ,jk)  )
312                        v_abl( ji, jj, jk, nt_a ) =                v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff
313                        zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_a ) - pgv_dta(ji  ,jj  ,jk)  )
314                        u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff )
315                        v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  v_abl( ji, jj, jk, nt_a )
316                     END DO
317                  END DO
318               END IF
319               !
320            !-------------
321            END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
322            !-------------
323
324         ENDIF ! ln_geos_winds
325
326      ENDIF ! SemiImp_Cor
327      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
328      !                            !  4 *** Advance u,v to time n+1
329      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
330      !
331      !  Vertical diffusion for u_abl
332      !-------------
333      DO jj = 1, jpj    ! outer loop
334      !-------------
335
336         DO jk = 3, jpkam1
337            DO ji = 1, jpi
338               z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal
339               z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal
340               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )  !       diagonal
341            END DO
342         END DO
343
344         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)
345            !++ Surface boundary condition
346            z_elem_a( ji, 2 ) = 0._wp
347            z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
348            !
349            zztmp1  = pcd_du(ji, jj)
350            zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) * rn_vfac
351#if defined key_si3
352            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
353            zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj) ) * rn_vfac
354            zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
355#endif
356            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
357            u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rDt_abl * zztmp2
358
359            ! idealized test cases only
360            !IF( ln_topbc_neumann ) THEN
361            !   !++ Top Neumann B.C.
362            !   z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
363            !   z_elem_c( ji,     jpka ) = 0._wp
364            !   z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
365            !   !u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl   ( ji, jj, jpka, nt_a )
366            !ELSE
367               !++ Top Dirichlet B.C.
368               z_elem_a( ji,     jpka )       = 0._wp
369               z_elem_c( ji,     jpka )       = 0._wp
370               z_elem_b( ji,     jpka )       = e3t_abl( jpka )
371               u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)
372            !ENDIF
373
374         END DO
375         !!
376         !! Matrix inversion
377         !! ----------------------------------------------------------
378         !DO ji = 2, jpi
379         DO ji = 1, jpi  !!GS: TBI
380            zcff                 =   1._wp / z_elem_b( ji, 2 )
381            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 )
382            u_abl (ji,jj,2,nt_a) =    zcff * u_abl(ji,jj,2,nt_a)
383         END DO
384
385         DO jk = 3, jpka
386            DO ji = 2, jpi
387               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
388               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
389               u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk  ,nt_a)   &
390               &          - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
391            END DO
392         END DO
393
394         DO jk = jpkam1,2,-1
395            DO ji = 2, jpi
396               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)
397            END DO
398         END DO
399
400      !-------------
401      END DO             ! end outer loop
402      !-------------
403
404      !
405      !  Vertical diffusion for v_abl
406      !-------------
407      DO jj = 2, jpj   ! outer loop
408      !-------------
409         !
410         DO jk = 3, jpkam1
411            DO ji = 1, jpi
412               z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
413               z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal
414               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal
415            END DO
416         END DO
417
418         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)
419            !++ Surface boundary condition
420            z_elem_a( ji, 2 ) = 0._wp
421            z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
422            !
423            zztmp1 = pcd_du(ji, jj)
424            zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) * rn_vfac
425#if defined key_si3
426            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
427            zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) ) * rn_vfac
428            zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
429#endif
430            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
431            v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
432
433            ! idealized test cases only
434            !IF( ln_topbc_neumann ) THEN
435            !   !++ Top Neumann B.C.
436            !   z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
437            !   z_elem_c( ji,     jpka ) = 0._wp
438            !   z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
439            !   !v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl   ( ji, jj, jpka, nt_a )
440            !ELSE
441               !++ Top Dirichlet B.C.
442               z_elem_a( ji,     jpka )       = 0._wp
443               z_elem_c( ji,     jpka )       = 0._wp
444               z_elem_b( ji,     jpka )       = e3t_abl( jpka )
445               v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)
446            !ENDIF
447
448         END DO
449         !!
450         !! Matrix inversion
451         !! ----------------------------------------------------------
452         DO ji = 1, jpi
453            zcff                 =  1._wp / z_elem_b( ji, 2 )
454            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       )
455            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a )
456         END DO
457
458         DO jk = 3, jpka
459            DO ji = 1, jpi
460               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
461               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
462               v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk  ,nt_a)   &
463               &          - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
464            END DO
465         END DO
466
467         DO jk = jpkam1,2,-1
468            DO ji = 1, jpi
469               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)
470            END DO
471         END DO
472         !
473      !-------------
474      END DO             ! end outer loop
475      !-------------
476
477      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
478      !                            !  5 *** Apply nudging on the dynamics and the tracers
479      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
480
481      IF( nn_dyn_restore > 0  ) THEN
482         !-------------
483         DO jk = 2, jpka    ! outer loop
484         !-------------
485            DO_2D_01_01
486               zcff1 = pblh( ji, jj )
487               zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
488               zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )
489               zmsk  = msk_abl(ji,jj)
490               zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   &
491                  &  + jp_alp1_dyn * zsig    + jp_alp0_dyn
492               zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points
493                                                             ! rn_Dt = rDt_abl / nn_fsbc
494               zcff  = zcff * rest_eq(ji,jj)
495               u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   &
496                  &                               + zcff   * pu_dta( ji, jj, jk       )
497               v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   &
498                  &                               + zcff   * pv_dta( ji, jj, jk       )
499            END_2D
500         !-------------
501         END DO             ! end outer loop
502         !-------------
503      END IF
504
505      !-------------
506      DO jk = 2, jpka    ! outer loop
507      !-------------
508         DO_2D_11_11
509            zcff1 = pblh( ji, jj )
510            zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
511            zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )
512            zmsk  = msk_abl(ji,jj)
513            zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   &
514               &  + jp_alp1_tra * zsig    + jp_alp0_tra
515            zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points
516                                                          ! rn_Dt = rDt_abl / nn_fsbc
517            tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   &
518               &                                       + zcff   * pt_dta( ji, jj, jk )
519
520            tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa )   &
521               &                                       + zcff   * pq_dta( ji, jj, jk )
522
523         END_2D
524      !-------------
525      END DO             ! end outer loop
526      !-------------
527      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
528      !                            !  6 *** MPI exchanges
529      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
530      !
531      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a)      , 'T', -1.                            )
532      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...
533      !
534#if defined key_iomput
535      ! 2D & first ABL level
536      IF ( iom_use("pblh"   ) ) CALL iom_put (    "pblh",    pblh(:,:             ) )
537      IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl",   u_abl(:,:,2,nt_a      ) )
538      IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl",   v_abl(:,:,2,nt_a      ) )
539      IF ( iom_use("tz1_abl") ) CALL iom_put ( "tz1_abl",  tq_abl(:,:,2,nt_a,jp_ta) )
540      IF ( iom_use("qz1_abl") ) CALL iom_put ( "qz1_abl",  tq_abl(:,:,2,nt_a,jp_qa) )
541      IF ( iom_use("uz1_dta") ) CALL iom_put ( "uz1_dta",  pu_dta(:,:,2           ) )
542      IF ( iom_use("vz1_dta") ) CALL iom_put ( "vz1_dta",  pv_dta(:,:,2           ) )
543      IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta",  pt_dta(:,:,2           ) )
544      IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta",  pq_dta(:,:,2           ) )
545      ! debug 2D
546      IF( ln_geos_winds ) THEN
547         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) )
548         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) )
549      END IF
550      IF( ln_hpgls_frc ) THEN
551         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo",  pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
552         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
553      END IF
554      ! 3D (all ABL levels)
555      IF ( iom_use("u_abl"   ) ) CALL iom_put ( "u_abl"   ,    u_abl(:,:,2:jpka,nt_a      ) )
556      IF ( iom_use("v_abl"   ) ) CALL iom_put ( "v_abl"   ,    v_abl(:,:,2:jpka,nt_a      ) )
557      IF ( iom_use("t_abl"   ) ) CALL iom_put ( "t_abl"   ,   tq_abl(:,:,2:jpka,nt_a,jp_ta) )
558      IF ( iom_use("q_abl"   ) ) CALL iom_put ( "q_abl"   ,   tq_abl(:,:,2:jpka,nt_a,jp_qa) )
559      IF ( iom_use("tke_abl" ) ) CALL iom_put ( "tke_abl" ,  tke_abl(:,:,2:jpka,nt_a      ) )
560      IF ( iom_use("avm_abl" ) ) CALL iom_put ( "avm_abl" ,  avm_abl(:,:,2:jpka           ) )
561      IF ( iom_use("avt_abl" ) ) CALL iom_put ( "avt_abl" ,  avt_abl(:,:,2:jpka           ) )
562      IF ( iom_use("mxlm_abl") ) CALL iom_put ( "mxlm_abl", mxlm_abl(:,:,2:jpka           ) )
563      IF ( iom_use("mxld_abl") ) CALL iom_put ( "mxld_abl", mxld_abl(:,:,2:jpka           ) )
564      ! debug 3D
565      IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta",  pu_dta(:,:,2:jpka) )
566      IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta",  pv_dta(:,:,2:jpka) )
567      IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta",  pt_dta(:,:,2:jpka) )
568      IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta",  pq_dta(:,:,2:jpka) )
569      IF( ln_geos_winds ) THEN
570         IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka) )
571         IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgv_dta(:,:,2:jpka) )
572      END IF
573      IF( ln_hpgls_frc ) THEN
574         IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo",  pgu_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) )
575         IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", -pgv_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) )
576      END IF
577#endif
578      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
579      !                            !  7 *** Finalize flux computation
580      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
581      !
582      DO_2D_11_11
583         ztemp          =  tq_abl( ji, jj, 2, nt_a, jp_ta )
584         zhumi          =  tq_abl( ji, jj, 2, nt_a, jp_qa )
585         zcff           = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
586         psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )   !GS: negative sign to respect aerobulk convention
587         pevp( ji, jj ) = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) )
588         rhoa( ji, jj ) = zcff
589      END_2D
590
591      DO_2D_01_01
592<<<<<<< .working
593         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * ( pssu(ji  ,jj) + pssu(ji-1,jj) ) 
594         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * ( pssv(ji,jj  ) + pssv(ji,jj-1) ) 
595||||||| .merge-left.r13136
596         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji  ,jj) + pssu(ji-1,jj) ) 
597         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj  ) + pssv(ji,jj-1) ) 
598=======
599         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji,jj) + pssu(ji-1,jj  ) )
600         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj) + pssv(ji  ,jj-1) )
601>>>>>>> .merge-right.r13211
602      END_2D
603      !
604      CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
605      !
606      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
607      DO_2D_11_11
608         zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   &
609            &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj) )   ! * msk_abl(ji,jj)
610         zztmp         = rhoa(ji,jj) * pcd_du(ji,jj)
611
612         pwndm (ji,jj) =         zcff
613         ptaum (ji,jj) = zztmp * zcff
614         zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
615         zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
616      END_2D
617      ! ... utau, vtau at U- and V_points, resp.
618      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
619      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
620      DO_2D_00_00
621         zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
622         zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
623         ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) )
624         zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
625         zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
626         ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) )
627      END_2D
628      !
629      CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1., ptauj(:,:), 'V', -1. )
630
631      CALL iom_put( "taum_oce", ptaum )
632
633      IF(sn_cfctl%l_prtctl) THEN
634         CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau   : ', mask1=umask,   &
635            &          tab2d_2=ptauj , clinfo2='          vtau   : ', mask2=vmask )
636         CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm   : ' )
637      ENDIF
638
639#if defined key_si3
640<<<<<<< .working
641      ! ------------------------------------------------------------ !
642      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
643      ! ------------------------------------------------------------ !
644      DO_2D_00_00           
645         ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   &
646            &                      * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
647         ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   &
648            &                      * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
649      END_2D
650      CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
651      !
652      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   &
653         &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' )
654||||||| .merge-left.r13136
655         ! ------------------------------------------------------------ !
656         !    Wind stress relative to the moving ice ( U10m - U_ice )   !
657         ! ------------------------------------------------------------ !
658         DO_2D_00_00
659           
660            zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) )
661            zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
662   
663            ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             &
664               &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          &
665               &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) )
666            ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             &
667               &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          &
668               &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) )
669         END_2D
670         CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
671         !
672         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   &
673            &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' )
674=======
675      ! ------------------------------------------------------------ !
676      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
677      ! ------------------------------------------------------------ !
678      DO_2D_00_00
679
680         zztmp1 = 0.5_wp * ( u_abl(ji+1,jj  ,2,nt_a) + u_abl(ji,jj,2,nt_a) )
681         zztmp2 = 0.5_wp * ( v_abl(ji  ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
682
683         ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             &
684            &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          &
685            &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) )
686         ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             &
687            &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          &
688            &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) )
689      END_2D
690      CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
691      !
692      IF(sn_cfctl%l_prtctl) THEN
693         CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask,   &
694            &          tab2d_2=ptauj_ice , clinfo2='          vtau_ice : ', mask2=vmask )
695      END IF
696>>>>>>> .merge-right.r13211
697#endif
698      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
699      !                            !  8 *** Swap time indices for the next timestep
700      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
701      nt_n = 1 + MOD( nt_n, 2)
702      nt_a = 1 + MOD( nt_a, 2)
703      !
704!---------------------------------------------------------------------------------------------------
705   END SUBROUTINE abl_stp
706!===================================================================================================
707
708
709
710
711!===================================================================================================
712   SUBROUTINE abl_zdf_tke( )
713!---------------------------------------------------------------------------------------------------
714
715      !!----------------------------------------------------------------------
716      !!                   ***  ROUTINE abl_zdf_tke  ***
717      !!
718      !! ** Purpose :   Time-step Turbulente Kinetic Energy (TKE) equation
719      !!
720      !! ** Method  : - source term due to shear
721      !!              - source term due to stratification
722      !!              - resolution of the TKE equation by inverting
723      !!                a tridiagonal linear system
724      !!
725      !! ** Action  : - en : now turbulent kinetic energy)
726      !!              - avmu, avmv : production of TKE by shear at u and v-points
727      !!                (= Kz dz[Ub] * dz[Un] )
728      !! ---------------------------------------------------------------------
729      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn
730      INTEGER, DIMENSION(1:jpi          )     ::   ikbl
731      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
732      REAL(wp)                                ::   zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2
733      REAL(wp)                                ::   zdU2, zdV2, zbuoy1, zbuoy2    ! zbuoy for BL89
734      REAL(wp)                                ::   zwndi, zwndj
735      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2
736      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2
737      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF
738      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a
739      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b
740      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_c
741      LOGICAL                                 ::   ln_Patankar    = .FALSE.
742      LOGICAL                                 ::   ln_dumpvar     = .FALSE.
743      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl
744      !
745      tind  = nt_n
746      ziRic = 1._wp / rn_Ric
747      ! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
748      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
749      !                            !  Advance TKE equation to time n+1
750      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
751      !-------------
752      DO jj = 1, jpj    ! outer loop
753      !-------------
754         !
755         ! Compute vertical shear
756         DO jk = 2, jpkam1
757            DO ji = 1, jpi
758               zcff        = 1.0_wp / e3w_abl( jk )**2
759               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
760               zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
761               zsh2(ji,jk) = zdU+zdV   !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 )
762            END DO
763         END DO
764         !
765         ! Compute brunt-vaisala frequency
766         DO jk = 2, jpkam1
767            DO ji = 1,jpi
768               zcff  = grav * itvref / e3w_abl( jk )
769               zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta)
770               zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        &
771                  &   - tq_abl( ji, jj, jk  , tind, jp_ta) * tq_abl( ji, jj, jk  , tind, jp_qa)
772               zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi)
773            END DO
774         END DO
775         !
776         ! Terms for the tridiagonal problem
777         DO jk = 2, jpkam1
778            DO ji = 1, jpi
779               zshear      = zsh2( ji, jk )                           ! zsh2 is already multiplied by Avm_abl at this point
780               zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
781               zbuoy       = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
782
783               z_elem_a( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk ) + Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk   ) ! lower-diagonal
784               z_elem_c( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk ) + Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal
785               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE
786                  z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
787                     &               + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)    ! diagonal
788                  tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )       ! right-hand-side
789               ELSE
790                  z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
791                     &               + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)   &  ! diagonal
792                     &               - e3w_abl(jk) * rDt_abl * zbuoy
793                  tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )                    ! right-hand-side
794               END IF
795            END DO
796         END DO
797
798         DO ji = 1,jpi    ! vector opt.
799            zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min )
800            zetop = tke_min
801
802            z_elem_a ( ji,     1       ) = 0._wp
803            z_elem_c ( ji,     1       ) = 0._wp
804            z_elem_b ( ji,     1       ) = 1._wp
805            tke_abl  ( ji, jj, 1, nt_a ) = zesrf
806
807            !++ Top Neumann B.C.
808            !z_elem_a ( ji,     jpka       ) = - 0.5 * rDt_abl * rn_Sch * (Avm_abl(ji,jj, jpka-1 )+Avm_abl(ji,jj, jpka ))  / e3t_abl( jpka )
809            !z_elem_c ( ji,     jpka       ) = 0._wp
810            !z_elem_b ( ji,     jpka       ) = e3w_abl(jpka) - z_elem_a(ji, jpka )
811            !tke_abl  ( ji, jj, jpka, nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n )
812
813            !++ Top Dirichlet B.C.
814            z_elem_a ( ji,     jpka       ) = 0._wp
815            z_elem_c ( ji,     jpka       ) = 0._wp
816            z_elem_b ( ji,     jpka       ) = 1._wp
817            tke_abl  ( ji, jj, jpka, nt_a ) = zetop
818
819            zbn2 ( ji, jj,    1 ) = zbn2 ( ji, jj,      2 )
820            zsh2 ( ji,        1 ) = zsh2 ( ji,          2 )
821            zbn2 ( ji, jj, jpka ) = zbn2 ( ji, jj, jpkam1 )
822            zsh2 ( ji,     jpka ) = zsh2 ( ji    , jpkam1 )
823         END DO
824         !!
825         !! Matrix inversion
826         !! ----------------------------------------------------------
827         DO ji = 1,jpi
828            zcff                  =  1._wp / z_elem_b( ji, 1 )
829            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       )
830            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a )
831         END DO
832
833         DO jk = 2, jpka
834            DO ji = 1,jpi
835               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
836               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
837               tke_abl(ji,jj,jk,nt_a) =   zcff * ( tke_abl(ji,jj,jk  ,nt_a)   &
838               &          - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
839            END DO
840         END DO
841
842         DO jk = jpkam1,1,-1
843            DO ji = 1,jpi
844               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)
845            END DO
846         END DO
847
848!!FL should not be needed because of Patankar procedure
849         tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
850
851         !!
852         !! Diagnose PBL height
853         !! ----------------------------------------------------------
854
855
856         !
857         ! arrays zRH, zFC and zCF are available at this point
858         ! and zFC(:, 1 ) = 0.
859         ! diagnose PBL height based on zsh2 and zbn2
860         zFC (  :  ,1) = 0._wp
861         ikbl( 1:jpi ) = 0
862
863         DO jk = 2,jpka
864            DO ji = 1, jpi
865               zcff  = ghw_abl( jk-1 )
866               zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
867               zcff  = ghw_abl( jk   )
868               zcff2 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
869               zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*(                 &
870                               zcff2 * ( zsh2( ji, jk  ) - ziRic * zbn2( ji, jj, jk   ) &
871                           - rn_Cek  * ( fft_abl( ji, jj  ) * fft_abl( ji, jj ) ) ) &
872                             + zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
873                           - rn_Cek  * ( fft_abl( ji, jj  ) * fft_abl( ji, jj ) ) ) &
874                           &                                                 )
875               IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
876            END DO
877         END DO
878         !
879         ! finalize the computation of the PBL height
880         DO ji = 1, jpi
881            jk = ikbl(ji)
882            IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
883               pblh( ji, jj ) =  (  ghw_abl( jk-1 ) * zFC( ji, jk   )       &
884                  &              -  ghw_abl( jk   ) * zFC( ji, jk-1 )       &
885                  &              ) / ( zFC( ji, jk   ) - zFC( ji, jk-1 ) )
886            ELSE IF( jk==2 ) THEN
887               pblh( ji, jj ) = ghw_abl(2   )
888            ELSE
889               pblh( ji, jj ) = ghw_abl(jpka)
890            END IF
891         END DO
892      !-------------
893      END DO
894      !-------------
895      !
896      ! Optional : could add pblh smoothing if pblh is noisy horizontally ...
897      IF(ln_smth_pblh) THEN
898         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing)
899         CALL smooth_pblh( pblh, msk_abl )
900         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing)
901      ENDIF
902      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
903      !                            !  Diagnostic mixing length computation
904      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
905      !
906      SELECT CASE ( nn_amxl )
907      !
908      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom
909#   define zlup zRH
910#   define zldw zFC
911         DO jj = 1, jpj     ! outer loop
912            !
913            DO ji = 1, jpi
914               mxld_abl( ji, jj,    1 ) = mxl_min
915               mxld_abl( ji, jj, jpka ) = mxl_min
916               mxlm_abl( ji, jj,    1 ) = mxl_min
917               mxlm_abl( ji, jj, jpka ) = mxl_min
918               zldw    ( ji,        1 ) = zrough(ji,jj) * rn_Lsfc
919               zlup    ( ji,     jpka ) = mxl_min
920            END DO
921            !
922            DO jk = 2, jpkam1
923               DO ji = 1, jpi
924                  zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
925                  mxlm_abl( ji, jj, jk ) = MAX( mxl_min,  &
926                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
927               END DO
928            END DO
929            !
930            ! Limit mxl
931            DO jk = jpkam1,1,-1
932               DO ji = 1, jpi
933                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
934               END DO
935            END DO
936            !
937            DO jk = 2, jpka
938               DO ji = 1, jpi
939                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
940               END DO
941            END DO
942            !
943!            DO jk = 1, jpka
944!               DO ji = 1, jpi
945!                  mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
946!                  mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ),  zlup( ji, jk ) )
947!               END DO
948!            END DO
949            !
950            DO jk = 1, jpka
951               DO ji = 1, jpi
952!                  zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
953                  zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
954                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
955                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
956               END DO
957            END DO
958            !
959         END DO
960#   undef zlup
961#   undef zldw
962         !
963         !
964      CASE ( 1 )           ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom
965#   define zlup zRH
966#   define zldw zFC
967         DO jj = 1, jpj     ! outer loop
968            !
969            DO jk = 2, jpkam1
970               DO ji = 1,jpi
971                              zcff        = 1.0_wp / e3w_abl( jk )**2
972                  zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
973                  zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
974                  zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
975                           ENDDO
976                        ENDDO
977                        !
978            DO ji = 1, jpi
979               zcff                      = zrough(ji,jj) * rn_Lsfc
980               mxld_abl ( ji, jj,    1 ) = zcff
981               mxld_abl ( ji, jj, jpka ) = mxl_min
982               mxlm_abl ( ji, jj,    1 ) = zcff
983               mxlm_abl ( ji, jj, jpka ) = mxl_min
984               zldw     ( ji,        1 ) = zcff
985               zlup     ( ji,     jpka ) = mxl_min
986            END DO
987            !
988            DO jk = 2, jpkam1
989               DO ji = 1, jpi
990                  zbuoy    = MAX( zbn2(ji, jj, jk), rsmall )
991                  zcff     = 2.*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) &
992                                &             + SQRT( rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.*zbuoy ) )
993                                  mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff )
994               END DO
995            END DO
996            !
997            ! Limit mxl
998            DO jk = jpkam1,1,-1
999               DO ji = 1, jpi
1000                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
1001               END DO
1002            END DO
1003            !
1004            DO jk = 2, jpka
1005               DO ji = 1, jpi
1006                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
1007               END DO
1008            END DO
1009            !
1010            DO jk = 1, jpka
1011               DO ji = 1, jpi
1012                  !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1013                  !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
1014                  zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1015                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
1016                  !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) )
1017                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
1018               END DO
1019            END DO
1020 !
1021         END DO
1022#   undef zlup
1023#   undef zldw
1024         !
1025      CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale
1026         !
1027#   define zlup zRH
1028#   define zldw zFC
1029! zCF is used for matrix inversion
1030!
1031       DO jj = 1, jpj      ! outer loop
1032
1033         DO ji = 1, jpi
1034            zcff             = zrough(ji,jj) * rn_Lsfc
1035            zlup( ji,    1 ) = zcff
1036            zldw( ji,    1 ) = zcff
1037            zlup( ji, jpka ) = mxl_min
1038            zldw( ji, jpka ) = mxl_min
1039         END DO
1040
1041         DO jk = 2,jpka-1
1042            DO ji = 1, jpi
1043               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
1044               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)
1045            END DO
1046         END DO
1047         !!
1048         !! BL89 search for lup
1049         !! ----------------------------------------------------------
1050         DO jk=2,jpka-1
1051            !
1052            DO ji = 1, jpi
1053               zCF(ji,1:jpka) = 0._wp
1054               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
1055               ln_foundl(ji ) = .false.
1056            END DO
1057            !
1058            DO jkup=jk+1,jpka-1
1059               DO ji = 1, jpi
1060                  zbuoy1 = MAX( zbn2(ji,jj,jkup  ), rsmall )
1061                  zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
1062                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
1063                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk)) &
1064                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )
1065                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1066                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
1067                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
1068                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
1069                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) )
1070                     zlup(ji,jk) = zcff
1071                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk)
1072                     ln_foundl(ji) = .true.
1073                  END IF
1074               END DO
1075            END DO
1076            !
1077         END DO
1078         !!
1079         !! BL89 search for ldwn
1080         !! ----------------------------------------------------------
1081         DO jk=2,jpka-1
1082            !
1083            DO ji = 1, jpi
1084               zCF(ji,1:jpka) = 0._wp
1085               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
1086               ln_foundl(ji ) = .false.
1087            END DO
1088            !
1089            DO jkdwn=jk-1,1,-1
1090               DO ji = 1, jpi
1091                  zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
1092                  zbuoy2 = MAX( zbn2(ji,jj,jkdwn  ), rsmall )
1093                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
1094                     &               * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
1095                                       + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) )
1096                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
1097                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
1098                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )
1099                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
1100                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) )
1101                     zldw(ji,jk) = zcff
1102                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  )
1103                     ln_foundl(ji) = .true.
1104                  END IF
1105               END DO
1106            END DO
1107            !
1108         END DO
1109
1110         DO jk = 1, jpka
1111            DO ji = 1, jpi
1112               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
1113               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1114               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
1115               mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
1116            END DO
1117         END DO
1118
1119      END DO
1120#   undef zlup
1121#   undef zldw
1122         !
1123     CASE ( 3 )           ! Bougeault & Lacarrere 89 length-scale
1124         !
1125#   define zlup zRH
1126#   define zldw zFC
1127! zCF is used for matrix inversion
1128!
1129       DO jj = 1, jpj      ! outer loop
1130          !
1131          DO jk = 2, jpkam1
1132             DO ji = 1,jpi
1133                            zcff        = 1.0_wp / e3w_abl( jk )**2
1134                zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
1135                zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
1136                zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
1137                         ENDDO
1138                  ENDDO
1139          zsh2(:,      1) = zsh2( :,      2)
1140          zsh2(:,   jpka) = zsh2( :, jpkam1)
1141
1142                 DO ji = 1, jpi
1143            zcff              = zrough(ji,jj) * rn_Lsfc
1144                        zlup( ji,    1 )  = zcff
1145            zldw( ji,    1 )  = zcff
1146            zlup( ji, jpka ) = mxl_min
1147            zldw( ji, jpka ) = mxl_min
1148         END DO
1149
1150         DO jk = 2,jpka-1
1151            DO ji = 1, jpi
1152               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
1153               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)
1154            END DO
1155         END DO
1156         !!
1157         !! BL89 search for lup
1158         !! ----------------------------------------------------------
1159         DO jk=2,jpka-1
1160            !
1161            DO ji = 1, jpi
1162               zCF(ji,1:jpka) = 0._wp
1163               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
1164               ln_foundl(ji ) = .false.
1165            END DO
1166            !
1167            DO jkup=jk+1,jpka-1
1168               DO ji = 1, jpi
1169                  zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall )
1170                  zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall )
1171                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) *                 &
1172                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk))      &
1173                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )    &
1174                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *        &
1175                                         &               ( SQRT(tke_abl( ji, jj, jkup  , nt_a ))*zsh2(ji,jkup  ) &
1176                                         &               + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) )
1177
1178                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
1179                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
1180                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
1181                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
1182                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) )
1183                     zlup(ji,jk) = zcff
1184                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk)
1185                     ln_foundl(ji) = .true.
1186                  END IF
1187               END DO
1188            END DO
1189            !
1190         END DO
1191         !!
1192         !! BL89 search for ldwn
1193         !! ----------------------------------------------------------
1194         DO jk=2,jpka-1
1195            !
1196            DO ji = 1, jpi
1197               zCF(ji,1:jpka) = 0._wp
1198               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
1199               ln_foundl(ji ) = .false.
1200            END DO
1201            !
1202            DO jkdwn=jk-1,1,-1
1203               DO ji = 1, jpi
1204                  zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
1205                  zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall )
1206                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
1207                     &               * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1))    &
1208                     &                 +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) )  &
1209                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *          &
1210                                         &               ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) &
1211                                         &               + SQRT(tke_abl( ji, jj, jkdwn  , nt_a ))*zsh2(ji,jkdwn  ) )
1212
1213                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
1214                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
1215                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )
1216                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
1217                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) )
1218                     zldw(ji,jk) = zcff
1219                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  )
1220                     ln_foundl(ji) = .true.
1221                  END IF
1222               END DO
1223            END DO
1224            !
1225         END DO
1226
1227         DO jk = 1, jpka
1228            DO ji = 1, jpi
1229               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
1230               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
1231               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
1232               mxld_abl( ji, jj, jk ) = MAX(  MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
1233            END DO
1234         END DO
1235
1236      END DO
1237#   undef zlup
1238#   undef zldw
1239         !
1240         !
1241      END SELECT
1242      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1243      !                            !  Finalize the computation of turbulent visc./diff.
1244      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1245
1246      !-------------
1247      DO jj = 1, jpj     ! outer loop
1248      !-------------
1249         DO jk = 1, jpka
1250            DO ji = 1, jpi  ! vector opt.
1251               zcff  = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk )  &
1252               &     * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) )
1253               zcff2 =  1. / ( 1. + zcff )   !<-- phi_z(z)
1254               zcff  = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
1255               !!FL: MAX function probably useless because of the definition of mxl_min
1256               Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   )
1257               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )
1258            END DO
1259         END DO
1260      !-------------
1261      END DO
1262      !-------------
1263
1264!---------------------------------------------------------------------------------------------------
1265   END SUBROUTINE abl_zdf_tke
1266!===================================================================================================
1267
1268
1269!===================================================================================================
1270   SUBROUTINE smooth_pblh( pvar2d, msk )
1271!---------------------------------------------------------------------------------------------------
1272
1273      !!----------------------------------------------------------------------
1274      !!                   ***  ROUTINE smooth_pblh  ***
1275      !!
1276      !! ** Purpose :   2D Hanning filter on atmospheric PBL height
1277      !!
1278      !! ---------------------------------------------------------------------
1279      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk
1280      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
1281      INTEGER                                     :: ji,jj
1282      REAL(wp)                                    :: smth_a, smth_b
1283      REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY
1284      REAL(wp)                                    :: zumsk,zvmsk
1285      !!
1286      !!=========================================================
1287      !!
1288      !! Hanning filter
1289      smth_a = 1._wp / 8._wp
1290      smth_b = 1._wp / 4._wp
1291      !
1292      DO_2D_11_10
1293         zumsk = msk(ji,jj) * msk(ji+1,jj)
1294         zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk
1295      END_2D
1296
1297      DO_2D_10_11
1298         zvmsk = msk(ji,jj) * msk(ji,jj+1)
1299         zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk
1300      END_2D
1301
1302      DO_2D_10_00
1303         zFY ( ji, jj  ) =   zdY ( ji, jj   )                        &
1304            & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   &
1305            &            -  (zdX ( ji, jj   ) - zdX( ji-1, jj   ))  )
1306      END_2D
1307
1308      DO_2D_00_10
1309         zFX( ji, jj  ) =    zdX( ji, jj   )                         &
1310           &    + smth_a*(  (zdY( ji+1, jj ) - zdY( ji+1, jj-1))     &
1311           &             -  (zdY( ji  , jj ) - zdY( ji  , jj-1)) )
1312      END_2D
1313
1314      DO_2D_00_00
1315         pvar2d( ji  ,jj ) = pvar2d( ji  ,jj )              &
1316  &         + msk(ji,jj) * smth_b * (                       &
1317  &                  zFX( ji, jj ) - zFX( ji-1, jj )        &
1318  &                 +zFY( ji, jj ) - zFY( ji, jj-1 )  )
1319      END_2D
1320
1321!---------------------------------------------------------------------------------------------------
1322   END SUBROUTINE smooth_pblh
1323!===================================================================================================
1324
1325!!======================================================================
1326END MODULE ablmod
Note: See TracBrowser for help on using the repository browser.