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

Last change on this file since 11322 was 11322, checked in by flemarie, 2 years ago

First implementation of ABL (see ticket #2131)

  • Update reference and cfg namelists for ORCA2_ICE_ABL
  • Run ABL over the ocean and BLK over sea-ice (ABL over sea-ice to come)
  • Bug fix in computation of pblh (+ add option to smooth pblh)
File size: 46.6 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      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)     ::   z_cft      !--FL--to be removed after the test phase   
87      !
88      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a
89      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_b
90      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_c
91      !
92      INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
93      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1
94      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv
95      !
96      !!---------------------------------------------------------------------     
97      !
98      IF(lwp .AND. kt == nit000) THEN                  ! control print
99         WRITE(numout,*)
100         WRITE(numout,*) 'abl_stp : ABL time stepping'
101         WRITE(numout,*) '~~~~~~'
102      ENDIF 
103      !
104      IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
105      !! Compute ustar squared as Cd || Uatm-Uoce ||^2
106      !! needed for surface boundary condition of TKE
107      !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk)
108      DO jj = 1,jpj
109         DO ji = 1,jpi
110            ustar2(ji,jj) = pCd_du(ji,jj) * pwndm(ji,jj)
111         END DO
112      END DO 
113      !
114      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
115      !                            !  1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
116      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
117
118      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
119         
120      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
121      !                            !  2 *** Advance tracers to time n+1
122      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
123     
124      !-------------
125      DO jj = 1, jpj    ! outer loop             !--> tq_abl computed on (1:jpi) x (1:jpj)
126      !-------------   
127         ! Compute matrix elements for interior points
128         DO jk = 3, jpkam1
129            DO ji = 1, jpi   ! vector opt.
130               z_elem_a( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )          ! lower-diagonal
131               z_elem_c( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )          ! upper-diagonal       
132               z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )  !       diagonal           
133            END DO
134         END DO 
135         ! Boundary conditions 
136         DO ji = 1, jpi   ! vector opt.
137            ! Neumann at the bottom         
138            z_elem_a( ji,     2              ) = 0._wp
139            z_elem_c( ji,     2              ) = - rdt * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                         
140            ! Homogeneous Neumann at the top
141            z_elem_a( ji,     jpka           ) = - rdt * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
142            z_elem_c( ji,     jpka           ) = 0._wp
143            z_elem_b( ji,     jpka           ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
144         END DO           
145                   
146         DO jtra = 1,jptq  ! loop on active tracers
147               
148            DO jk = 3, jpkam1
149               DO ji = fs_2, fs_jpim1                 
150                  tq_abl  ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl  ( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side
151               END DO
152            END DO
153
154            IF(jtra == jp_ta) THEN
155               DO ji = 1,jpi  ! surface boundary condition for temperature             
156                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * psen(ji, jj) 
157                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   &
158                    &                               + rdt * psen(ji, jj) * ( psst(ji, jj) + rt0 )               
159                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
160               END DO
161            ELSE   
162               DO ji = 1,jpi  ! surface boundary condition for humidity             
163                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * pevp(ji, jj) 
164                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   &
165                    &                               + rdt * pevp(ji, jj) * pssq(ji, jj)               
166                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
167               END DO                                     
168            END IF
169            !!
170            !! Matrix inversion
171            !! ----------------------------------------------------------
172            DO ji = 1,jpi           
173               zcff                       =  1._wp / z_elem_b( ji, 2 )
174               zCF   (ji,   2           ) = - zcff * z_elem_c( ji, 2 )
175               tq_abl(ji,jj,2,nt_a,jtra) =    zcff * tq_abl(ji,jj,2,nt_a,jtra)
176            END DO             
177
178            DO jk = 3, jpka           
179               DO ji = 1,jpi           
180                  zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
181                  zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
182                  tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk  ,nt_a,jtra)   &
183                  &                - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) )
184               END DO
185            END DO
186!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...         
187            DO jk = jpkam1,2,-1           
188               DO ji = 1,jpi 
189                  tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) +    &
190                     &                        zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
191               END DO
192            END DO
193         
194         END DO   !<-- loop on tracers                   
195         !!
196      !-------------
197      END DO             ! end outer loop
198      !-------------   
199     
200     
201      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
202      !                            !  3 *** Compute Coriolis term with geostrophic guide
203      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
204      !-------------       
205      DO jk = 2, jpka    ! outer loop
206      !-------------
207         !             
208         ! Advance u_abl & v_abl to time n+1
209         DO jj = 1, jpj
210            DO ji = 1, jpi           
211               zcff = ( fft_abl(ji,jj) * rdt )*( fft_abl(ji,jj) * rdt )  ! (f dt)**2
212     
213               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  &
214                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    &
215                  &                 +  rdt * fft_abl(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  &
216                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff)
217                 
218               v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  &
219                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   &
220                  &                 -  rdt * fft_abl(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) &
221                  &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)               
222            END DO
223         END DO 
224         !                                   
225      !-------------
226      END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
227      !-------------               
228      !
229      IF( ln_geos_winds ) THEN
230         DO jj = 1, jpj    ! outer loop       
231            DO jk = 1, jpka
232               DO ji = 1, jpi 
233                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
234                     &                      - rdt * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
235                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   &
236                     &                      + rdt * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
237               END DO
238            END DO
239         END DO     
240      END IF 
241      !-------------           
242      !
243      IF( ln_hpgls_frc ) THEN
244         DO jj = 1, jpj    ! outer loop       
245            DO jk = 1, jpka
246               DO ji = 1, jpi 
247                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgu_dta(ji,jj,jk)
248                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgv_dta(ji,jj,jk) 
249               ENDDO
250            ENDDO
251         ENDDO 
252      END IF
253
254      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
255      !                            !  4 *** Advance u,v to time n+1
256      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
257      !
258      !  Vertical diffusion for u_abl     
259      !-------------
260      DO jj = 1, jpj    ! outer loop
261      !-------------   
262
263         DO jk = 3, jpkam1
264            DO ji = 1, jpi 
265               z_elem_a( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal
266               z_elem_c( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal               
267               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal
268            END DO
269         END DO 
270           
271         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)           
272            !++ Surface boundary condition
273            z_elem_a( ji,     2    ) = 0._wp
274            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
275            zcff  = rdt * pcd_du(ji, jj)
276            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff 
277            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) ) 
278            !++ Top Neumann B.C.
279            !z_elem_a( ji,     jpka ) = - 0.5_wp * rdt * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )
280            !z_elem_c( ji,     jpka ) = 0._wp
281            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
282            !++ Top Dirichlet B.C.
283            z_elem_a( ji,     jpka )       = 0._wp
284            z_elem_c( ji,     jpka )       = 0._wp
285            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
286            u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)             
287         END DO 
288         !!
289         !! Matrix inversion
290         !! ----------------------------------------------------------
291         DO ji = 2, jpi         
292            zcff                 =   1._wp / z_elem_b( ji, 2 )
293            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 ) 
294            u_abl (ji,jj,2,nt_a) =    zcff * u_abl(ji,jj,2,nt_a)
295         END DO             
296
297         DO jk = 3, jpka           
298            DO ji = 2, jpi             
299               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
300               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
301               u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk  ,nt_a)   &
302               &          - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
303            END DO
304         END DO
305           
306         DO jk = jpkam1,2,-1           
307            DO ji = 2, jpi
308               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)
309            END DO
310         END DO
311                                         
312      !-------------
313      END DO             ! end outer loop
314      !-------------   
315
316      !
317      !  Vertical diffusion for v_abl     
318      !-------------
319      DO jj = 2, jpj   ! outer loop
320      !-------------   
321         !
322         DO jk = 3, jpkam1
323            DO ji = 1, jpi   
324               z_elem_a( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
325               z_elem_c( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal             
326               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal
327            END DO
328         END DO
329
330         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)         
331            !++ Surface boundary condition
332            z_elem_a( ji,     2    ) = 0._wp
333            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
334            zcff  = rdt * pcd_du(ji, jj)
335            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff 
336            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) ) 
337            !++ Top Neumann B.C.           
338            !z_elem_a( ji,     jpka ) = -rdt * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
339            !z_elem_c( ji,     jpka ) = 0._wp
340            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
341            !++ Top Dirichlet B.C.
342            z_elem_a( ji,     jpka )       = 0._wp
343            z_elem_c( ji,     jpka )       = 0._wp
344            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
345            v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)             
346         END DO 
347         !!
348         !! Matrix inversion
349         !! ----------------------------------------------------------
350         DO ji = 1, jpi             
351            zcff                 =  1._wp / z_elem_b( ji, 2 )
352            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       ) 
353            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a ) 
354         END DO             
355
356         DO jk = 3, jpka           
357            DO ji = 1, jpi                   
358               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
359               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
360               v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk  ,nt_a)   &
361               &          - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
362            END DO
363         END DO
364           
365         DO jk = jpkam1,2,-1           
366            DO ji = 1, jpi   
367               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)
368            END DO
369         END DO
370         !                                         
371      !-------------
372      END DO             ! end outer loop
373      !-------------   
374 
375      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
376      !                            !  5 *** Apply nudging on the dynamics and the tracers
377      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
378      z_cft(:,:,:) = 0._wp     
379     
380      IF( nn_dyn_restore > 0  ) THEN
381         !-------------
382         DO jk = 2, jpka    ! outer loop
383         !-------------       
384            DO jj = 2, jpj
385               DO ji = 2, jpi
386                  zcff1  = pblh( ji, jj )
387                  zsig   = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )                       
388                  zsig  =            MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
389                  zmsk  = msk_abl(ji,jj)
390                  zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   &
391                     &  + jp_alp1_dyn * zsig    + jp_alp0_dyn
392                  zcff  = (1._wp-zmsk) + zmsk * rdt * zcff2   ! zcff = 1 for masked points
393                 
394              zcff  = zcff * rest_eq(ji,jj) ; z_cft( ji, jj, jk ) = zcff
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      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
508      !                            !  8 *** Swap time indices for the next timestep
509      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
510      nt_n = 1 + MOD( kt  , 2)
511      nt_a = 1 + MOD( kt+1, 2)
512      !   
513!---------------------------------------------------------------------------------------------------
514   END SUBROUTINE abl_stp
515!===================================================================================================
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533!===================================================================================================
534   SUBROUTINE abl_zdf_tke( )
535!---------------------------------------------------------------------------------------------------
536
537      !!----------------------------------------------------------------------
538      !!                   ***  ROUTINE abl_zdf_tke  ***
539      !!
540      !! ** Purpose :   Time-step Turbulente Kinetic Energy (TKE) equation
541      !!
542      !! ** Method  : - source term due to shear
543      !!              - source term due to stratification
544      !!              - resolution of the TKE equation by inverting
545      !!                a tridiagonal linear system
546      !!
547      !! ** Action  : - en : now turbulent kinetic energy)
548      !!              - avmu, avmv : production of TKE by shear at u and v-points
549      !!                (= Kz dz[Ub] * dz[Un] )
550      !! ---------------------------------------------------------------------
551      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn 
552      INTEGER, DIMENSION(1:jpi          )     ::   ikbl
553      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
554      REAL(wp)                                ::   zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2
555      REAL(wp)                                ::   zdU2,zdV2     
556      REAL(wp)                                ::   zwndi,zwndj
557      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2
558      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2
559      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF       
560      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a
561      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b
562      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_c
563      LOGICAL                                 ::   ln_Patankar    = .FALSE.
564      LOGICAL                                 ::   ln_dumpvar     = .FALSE.
565      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl 
566      !
567      tind  = nt_n
568      ziRic = 1._wp / rn_Ric
569      ! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
570      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
571      !                            !  Advance TKE equation to time n+1
572      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
573      !-------------
574      DO jj = 1, jpj    ! outer loop
575      !-------------
576         !
577         ! Compute vertical shear         
578         DO jk = 2, jpkam1
579            DO ji = 1,jpi 
580               zcff        = 1.0_wp / e3w_abl( jk )**2               
581               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2                           
582               zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
583               zsh2(ji,jk) = zdU+zdV
584            END DO
585         END DO   
586         !
587         ! Compute brunt-vaisala frequency
588         DO jk = 2, jpkam1
589            DO ji = 1,jpi 
590               zcff  = grav * itvref / e3w_abl( jk ) 
591               zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta)
592               zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        &
593                  &   - tq_abl( ji, jj, jk  , tind, jp_ta) * tq_abl( ji, jj, jk  , tind, jp_qa)
594               zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi)
595            END DO
596         END DO 
597         !
598         ! Terms for the tridiagonal problem
599         DO jk = 2, jpkam1
600            DO ji = 1,jpi 
601               zshear       =                 zsh2( ji,     jk )   ! zsh2 is already multiplied by Avm_abl at this point
602               zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation         
603               zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
604               
605               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
606               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         
607               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE
608                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
609                     &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal       
610                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * ( zbuoy + zshear ) )             ! right-hand-side
611               ELSE
612                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
613                     &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal   
614                     &                     - e3w_abl(jk) * rdt * zbuoy   
615                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt *  zshear )             ! right-hand-side                     
616               END IF
617            END DO
618         END DO 
619                           
620         DO ji = 1,jpi    ! vector opt.           
621            zesrf   =  MAX( 4.63_wp * ustar2(ji,jj), tke_min )           
622            zetop   = tke_min             
623            z_elem_a ( ji,     1    ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp                                                   
624            z_elem_a ( ji,     jpka ) = 0._wp; z_elem_c ( ji,     jpka ) = 0._wp; z_elem_b ( ji,     jpka ) = 1._wp           
625            tke_abl( ji, jj,    1, nt_a ) = zesrf
626            tke_abl( ji, jj, jpka, nt_a ) = zetop 
627            zbn2(ji,jj,   1) = zbn2( ji,jj,   2) 
628            zsh2(ji,      1) = zsh2( ji,      2)                                     
629            zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 
630            zsh2(ji,   jpka) = zsh2( ji  , jpkam1)
631         END DO       
632         !!
633         !! Matrix inversion
634         !! ----------------------------------------------------------
635         DO ji = 1,jpi
636            zcff                  =  1._wp / z_elem_b( ji, 1 )
637            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       ) 
638            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a ) 
639         END DO             
640
641         DO jk = 2, jpka           
642            DO ji = 1,jpi
643               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
644               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
645               tke_abl(ji,jj,jk,nt_a) =   zcff * ( tke_abl(ji,jj,jk  ,nt_a)   &
646               &          - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
647            END DO
648         END DO
649           
650         DO jk = jpkam1,1,-1           
651            DO ji = 1,jpi
652               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)
653            END DO
654         END DO
655         
656!!FL should not be needed because of Patankar procedure 
657         tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
658
659         !!
660         !! Diagnose PBL height
661         !! ----------------------------------------------------------
662         
663         
664         !                                                       
665         ! arrays zRH, zFC and zCF are available at this point
666         ! and zFC(:, 1 ) = 0.
667         ! diagnose PBL height based on zsh2 and zbn2
668         zFC (  :  ,1) = 0._wp
669         ikbl( 1:jpi ) = 0 
670         
671         DO jk = 2,jpka
672            DO ji = 1, jpi 
673               zcff  = ghw_abl( jk-1 )
674               zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
675               zcff  = ghw_abl( jk   )
676               zcff2 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
677               zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*(                 &
678                               zcff2 * ( zsh2( ji, jk  ) - ziRic * zbn2( ji, jj, jk   ) &
679                           - rn_Cek  * ( fft_abl( ji, jj  ) * fft_abl( ji, jj ) ) ) &
680                             + zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
681                           - rn_Cek  * ( fft_abl( ji, jj  ) * fft_abl( ji, jj ) ) ) &
682                           &                                                 )
683               IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
684            END DO
685         END DO
686         !
687         ! finalize the computation of the PBL height
688         DO ji = 1, jpi
689            jk = ikbl(ji)
690            IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
691               pblh( ji, jj ) =  (  ghw_abl( jk-1 ) * zFC( ji, jk   )       &
692                  &              -  ghw_abl( jk   ) * zFC( ji, jk-1 )       &
693                  &              ) / ( zFC( ji, jk   ) - zFC( ji, jk-1 ) )
694            ELSE IF( jk==2 ) THEN
695               pblh( ji, jj ) = ghw_abl(2   )
696            ELSE
697               pblh( ji, jj ) = ghw_abl(jpka)
698            END IF
699         END DO
700         ! Optional : could add pblh smoothing if pblh is noisy horizontally ...
701         !
702      !-------------     
703      END DO     
704      !-------------
705     IF(ln_smth_pblh) THEN
706        CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
707        CALL smooth_pblh( pblh, msk_abl )
708        CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)   
709     ENDIF
710      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
711      !                            !  Diagnostic mixing length computation
712      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
713      !
714      SELECT CASE ( nn_amxl )
715      !
716      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom         
717#   define zlup zRH     
718#   define zldw zFC 
719         DO jj = 1, jpj     ! outer loop
720            !
721            DO ji = 1, jpi
722               mxl_abl ( ji, jj,    1 )  = 0._wp
723               mxl_abl ( ji, jj, jpka )  = mxl_min
724               zldw( ji,        1 )  = 0._wp             
725               zlup( ji,     jpka )  = 0._wp
726            END DO
727            !
728            DO jk = 2, jpkam1   
729               DO ji = 1, jpi 
730                  zbuoy             = MAX( zbn2(ji, jj, jk), rsmall ) 
731                  mxl_abl( ji, jj, jk ) = MAX( mxl_min,  &
732                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   ) 
733               END DO
734            END DO   
735            !
736            ! Limit mxl
737            DO jk = jpkam1,1,-1   
738               DO ji = 1, jpi 
739                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) )
740               END DO
741            END DO   
742            !
743            DO jk = 2, jpka
744               DO ji = 1, jpi 
745                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )   
746               END DO
747            END DO 
748            !
749            DO jk = 1, jpka
750               DO ji = 1, jpi
751                  mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )   
752               END DO
753            END DO                         
754            !
755         END DO
756#   undef zlup     
757#   undef zldw 
758         !
759         !
760      CASE ( 1 )             ! length-scale computed as the distance to the PBL height
761         DO jj = 1,jpj      ! outer loop
762            !
763            DO ji = 1, jpi   ! vector opt.
764               zcff = 1._wp / pblh( ji, jj )     ! inverse of hbl             
765               DO jk = 1, jpka             
766                  zsig  = MIN( zcff * ghw_abl( jk ), 1. )   
767                  zcff1 = pblh( ji, jj )                 
768                  mxl_abl( ji, jj, jk ) =  mxl_min                           &
769                     &  +  zsig         * (  amx1*zcff1 + bmx1*mxl_min ) &
770                     &  +  zsig *  zsig * (  amx2*zcff1 + bmx2*mxl_min ) &
771                     &  +  zsig**3      * (  amx3*zcff1 + bmx3*mxl_min ) &
772                     &  +  zsig**4      * (  amx4*zcff1 + bmx4*mxl_min ) &
773                     &  +  zsig**5      * (  amx5*zcff1 + bmx5*mxl_min )
774               END DO
775            END DO 
776            !           
777         END DO           
778         !
779      CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale
780         !
781#   define zlup zRH     
782#   define zldw zFC           
783! zCF is used for matrix inversion
784!             
785       DO jj = 1, jpj      ! outer loop
786         
787         DO ji = 1, jpi           
788            zlup( ji,    1 ) = mxl_min
789            zldw( ji,    1 ) = mxl_min               
790            zlup( ji, jpka ) = mxl_min
791            zldw( ji, jpka ) = mxl_min               
792         END DO           
793         
794         DO jk = 2,jpka-1
795            DO ji = 1, jpi
796               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
797               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)       
798            END DO
799         END DO         
800         !!
801         !! BL89 search for lup
802         !! ----------------------------------------------------------           
803         DO jk=2,jpka-1 
804            !
805            DO ji = 1, jpi
806               zCF(ji,1:jpka) = 0._wp
807               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
808               ln_foundl(ji ) = .false.
809            END DO   
810            !           
811            DO jkup=jk+1,jpka-1
812               DO ji = 1, jpi
813                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
814                     &               ( zbn2(ji,jj,jkup  )*(ghw_abl(jkup  )-ghw_abl(jk)) &
815                     &               + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )
816                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
817                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
818                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)                   
819                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
820                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
821                     zlup(ji,jk) = zcff               
822                     ln_foundl(ji) = .true.
823                  END IF
824               END DO
825            END DO
826            !
827         END DO   
828         !!
829         !! BL89 search for ldwn
830         !! ----------------------------------------------------------         
831         DO jk=2,jpka-1         
832            !
833            DO ji = 1, jpi
834               zCF(ji,1:jpka) = 0._wp
835               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
836               ln_foundl(ji ) = .false.
837            END DO 
838            !   
839            DO jkdwn=jk-1,1,-1
840               DO ji = 1, jpi             
841                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
842                     &               * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
843                                       + zbn2(ji,jj,jkdwn  )*(ghw_abl(jk)-ghw_abl(jkdwn  )) )   
844                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
845                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
846                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )             
847                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
848                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
849                     zldw(ji,jk) = zcff 
850                     ln_foundl(ji) = .true.               
851                  END IF                                   
852               END DO
853            END DO
854            !     
855         END DO
856
857         DO jk = 1, jpka
858            DO ji = 1, jpi 
859               mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min ) 
860            END DO
861         END DO 
862
863      END DO
864#   undef zlup     
865#   undef zldw           
866         !
867      END SELECT     
868      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
869      !                            !  Finalize the computation of turbulent visc./diff.
870      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
871     
872      !-------------
873      DO jj = 1, jpj     ! outer loop
874      !-------------
875         DO jk = 1, jpka   
876            DO ji = 1, jpi  ! vector opt.
877               zcff              = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk )  &
878                  &                                       * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) ) 
879               zcff2             =  1. / ( 1. + zcff )   !<-- phi_z(z)
880               zcff              = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
881!!FL: MAX function probably useless because of the definition of mxl_min             
882               Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   )
883               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )                             
884            END DO
885         END DO
886      !-------------         
887      END DO     
888      !-------------
889
890!---------------------------------------------------------------------------------------------------
891   END SUBROUTINE abl_zdf_tke
892!===================================================================================================
893
894
895!===================================================================================================
896   SUBROUTINE smooth_pblh( pvar2d, msk )
897!---------------------------------------------------------------------------------------------------
898
899      !!----------------------------------------------------------------------
900      !!                   ***  ROUTINE smooth_pblh  ***
901      !!
902      !! ** Purpose :   2D Hanning filter on atmospheric PBL height
903      !!
904      !! ---------------------------------------------------------------------
905     REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk   
906     REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
907      INTEGER                                     :: ji,jj
908     REAL(wp)                                    :: smth_a, smth_b
909     REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY
910     REAL(wp)                                    :: zumsk,zvmsk
911      !!
912      !!=========================================================
913      !!
914      !! Hanning filter
915      smth_a = 1._wp / 8._wp
916      smth_b = 1._wp / 4._wp
917      !
918      DO jj=1,jpj
919         DO ji=1,jpi-1
920            zumsk = msk(ji,jj) * msk(ji+1,jj)
921            zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk
922         END DO
923      END DO     
924     
925     DO jj=1,jpj-1
926         DO ji=1,jpi
927            zvmsk = msk(ji,jj) * msk(ji,jj+1)
928            zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk
929         END DO
930      END DO
931     
932     DO jj=1,jpj-1
933         DO ji=2,jpi-1
934            zFY ( ji, jj  ) =   zdY ( ji, jj   )                        &
935               & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   &
936               &            -  (zdX ( ji, jj   ) - zdX( ji-1, jj   ))  )
937         END DO
938      END DO
939
940      DO jj=2,jpj-1
941         DO ji=1,jpi-1
942            zFX( ji, jj  ) =    zdX( ji, jj   )                         &
943              &    + smth_a*(  (zdY( ji+1, jj ) - zdY( ji+1, jj-1))     &
944              &             -  (zdY( ji  , jj ) - zdY( ji  , jj-1)) )
945         END DO
946      END DO     
947
948     DO jj = 2, jpj-1
949         DO ji = 2,jpi-1
950            pvar2d( ji  ,jj ) = pvar2d( ji  ,jj )              &
951     &         + msk(ji,jj) * smth_b * (                       &
952     &                  zFX( ji, jj ) - zFX( ji-1, jj )        &
953     &                 +zFY( ji, jj ) - zFY( ji, jj-1 )  )
954         END DO
955      END DO 
956     !!
957!---------------------------------------------------------------------------------------------------
958   END SUBROUTINE smooth_pblh
959!===================================================================================================
960
961!!======================================================================
962END MODULE ablmod
Note: See TracBrowser for help on using the repository browser.