source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ABL/ablmod.F90 @ 12468

Last change on this file since 12468 was 12468, checked in by davestorkey, 9 months ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2:

  1. Alter ABL, ICE and TOP timestepping variables to be consistent with new schema: rdt_abl —> rDt_abl rdt_ice —> rDt_ice r1_rdt_ice —> r1_Dt_ice rdttrc —> rn_Dt (always equal to ocean timestep parameter in namelist) r2dttrc —> rDt_trc (current tracer timestep)
  2. Reinstate rn_scal_load (revert previous change): rn_load —> rn_scal_load

Passes SETTE and bit compares with the trunk@12436.

File size: 52.2 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_?fac
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
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,pssq_ice,pcd_du_ice  & 
47              &          , psen_ice, pevp_ice, pwndm_ice, pfrac_oce      & 
48              &          , ptaui_ice, ptauj_ice                          &
49#endif           
50           &   )
51!---------------------------------------------------------------------------------------------------
52
53      !!---------------------------------------------------------------------
54      !!                    ***  ROUTINE abl_stp ***
55      !!
56      !! ** Purpose :   Time-integration of the ABL model
57      !!
58      !! ** Method  :   Compute atmospheric variables : vertical turbulence
59      !!                             + Coriolis term + newtonian relaxation
60      !!               
61      !! ** Action  : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
62      !!              - Advance tracers to time n+1 (Euler backward scheme)
63      !!              - Compute Coriolis term with forward-backward scheme (possibly with geostrophic guide)
64      !!              - Advance u,v to time n+1 (Euler backward scheme)
65      !!              - Apply newtonian relaxation on the dynamics and the tracers
66      !!              - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
67      !!
68      !!----------------------------------------------------------------------
69      INTEGER  , INTENT(in   )                   ::   kt         ! time step index
70      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psst       ! sea-surface temperature [Celsius]
71      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu       ! sea-surface u (U-point)
72      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point)
73      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq       ! sea-surface humidity
74      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pu_dta     ! large-scale windi
75      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pv_dta     ! large-scale windj
76      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgu_dta    ! large-scale hpgi
77      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgv_dta    ! large-scale hpgj
78      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pt_dta     ! large-scale pot. temp.
79      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pq_dta     ! large-scale humidity
80      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp_dta   ! sea-level pressure
81      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du     ! Cd x Du (T-point)
82      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   psen       ! Ch x Du
83      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pevp       ! Ce x Du
84      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd||     
85      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui      ! taux
86      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj      ! tauy
87      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau||
88      !
89#if defined key_si3   
90      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptm_su       ! ice-surface temperature [K]
91      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu_ice     ! ice-surface u (U-point)
92      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point)     
93      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity
94      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du_ice   ! Cd x Du over ice (T-point)
95      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psen_ice     ! Ch x Du over ice (T-point)
96      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pevp_ice     ! Ce x Du over ice (T-point)
97      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice|| 
98      !REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pfrac_oce     !!GS: out useless ?
99      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pfrac_oce      !
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      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)    ::   z_cft      !--FL--to be removed after the test phase   
107      !
108      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a
109      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_b
110      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_c
111      !
112      INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
113      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
114      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
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      !! Compute ustar squared as Cd || Uatm-Uoce ||^2
126      !! needed for surface boundary condition of TKE
127      !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk)
128      DO_2D_11_11
129         zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj)
130#if defined key_si3
131         zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 
132         ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 
133#else
134         ustar2(ji,jj) = zzoce   
135#endif
136      END_2D
137      !
138      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
139      !                            !  1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
140      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
141
142      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
143         
144      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
145      !                            !  2 *** Advance tracers to time n+1
146      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
147     
148      !-------------
149      DO jj = 1, jpj    ! outer loop             !--> tq_abl computed on (1:jpi) x (1:jpj)
150      !-------------   
151         ! Compute matrix elements for interior points
152         DO jk = 3, jpkam1
153            DO ji = 1, jpi   ! vector opt.
154               z_elem_a( ji,     jk              ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
155               z_elem_c( ji,     jk              ) = - rDt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal       
156               z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal           
157            END DO
158         END DO 
159         ! Boundary conditions 
160         DO ji = 1, jpi   ! vector opt.
161            ! Neumann at the bottom         
162            z_elem_a( ji,     2              ) = 0._wp
163            z_elem_c( ji,     2              ) = - rDt_abl * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                         
164            ! Homogeneous Neumann at the top
165            z_elem_a( ji,     jpka           ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
166            z_elem_c( ji,     jpka           ) = 0._wp
167            z_elem_b( ji,     jpka           ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
168         END DO           
169
170         DO jtra = 1,jptq  ! loop on active tracers
171               
172            DO jk = 3, jpkam1
173               DO ji = 1,jpi
174                  tq_abl  ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl  ( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side
175               END DO
176            END DO
177
178            IF(jtra == jp_ta) THEN
179               DO ji = 1,jpi  ! boundary conditions for temperature             
180                  zztmp1 = psen(ji, jj) 
181                  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 ) 
182#if defined key_si3             
183                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 
184                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
185#endif             
186                  z_elem_b( ji,     2                ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1             
187                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rDt_abl * zztmp2               
188                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
189               END DO
190            ELSE   
191               DO ji = 1,jpi  ! boundary conditions for humidity             
192                  zztmp1 = pevp(ji, jj) 
193                  zztmp2 = pevp(ji, jj) * pssq(ji, jj)   
194#if defined key_si3             
195                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 
196                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)   
197#endif   
198                  z_elem_b( ji,     2                ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
199                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rDt_abl * zztmp2               
200                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
201               END DO                                     
202            END IF
203            !!
204            !! Matrix inversion
205            !! ----------------------------------------------------------
206            DO ji = 1,jpi           
207               zcff                       =  1._wp / z_elem_b( ji, 2 )
208               zCF   (ji,   2           ) = - zcff * z_elem_c( ji, 2 )
209               tq_abl(ji,jj,2,nt_a,jtra) =    zcff * tq_abl(ji,jj,2,nt_a,jtra)
210            END DO             
211
212            DO jk = 3, jpka           
213               DO ji = 1,jpi           
214                  zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
215                  zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
216                  tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk  ,nt_a,jtra)   &
217                  &                - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) )
218               END DO
219            END DO
220!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...         
221            DO jk = jpkam1,2,-1           
222               DO ji = 1,jpi 
223                  tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) +    &
224                     &                        zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
225               END DO
226            END DO
227         
228         END DO   !<-- loop on tracers                   
229         !!
230      !-------------
231      END DO             ! end outer loop
232      !-------------   
233     
234     
235      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
236      !                            !  3 *** Compute Coriolis term with geostrophic guide
237      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
238      !-------------       
239      DO jk = 2, jpka    ! outer loop
240      !-------------
241         !             
242         ! Advance u_abl & v_abl to time n+1
243         DO_2D_11_11
244            zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2
245   
246            u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  &
247               &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    &
248               &                 +  rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  &
249               &                               / (1._wp + gamma_Cor*gamma_Cor*zcff)
250               
251            v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  &
252               &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   &
253               &                 -  rDt_abl * fft_abl(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) &
254               &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)               
255         END_2D
256         !                                   
257      !-------------
258      END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
259      !-------------               
260      !
261      IF( ln_geos_winds ) THEN
262         DO jj = 1, jpj    ! outer loop       
263            DO jk = 1, jpka
264               DO ji = 1, jpi 
265                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
266                     &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
267                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   &
268                     &                      + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
269               END DO
270            END DO
271         END DO     
272      END IF 
273      !-------------           
274      !
275      IF( ln_hpgls_frc ) THEN
276         DO jj = 1, jpj    ! outer loop       
277            DO jk = 1, jpka
278               DO ji = 1, jpi 
279                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk)
280                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 
281               ENDDO
282            ENDDO
283         ENDDO 
284      END IF
285     
286      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
287      !                            !  4 *** Advance u,v to time n+1
288      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
289      !
290      !  Vertical diffusion for u_abl     
291      !-------------
292      DO jj = 1, jpj    ! outer loop
293      !-------------   
294
295         DO jk = 3, jpkam1
296            DO ji = 1, jpi 
297               z_elem_a( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal
298               z_elem_c( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal               
299               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal
300            END DO
301         END DO 
302           
303         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)           
304            !++ Surface boundary condition
305            z_elem_a( ji,     2    ) = 0._wp
306            z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
307            !
308         zztmp1  = pcd_du(ji, jj)
309         zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) )       
310#if defined key_si3             
311            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
312         zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) )
313         zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
314#endif           
315         z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1         
316            u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rDt_abl * zztmp2
317         
318            !++ Top Neumann B.C.
319            !z_elem_a( ji,     jpka ) = - 0.5_wp * rDt_abl * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )
320            !z_elem_c( ji,     jpka ) = 0._wp
321            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
322            !++ Top Dirichlet B.C.
323            z_elem_a( ji,     jpka )       = 0._wp
324            z_elem_c( ji,     jpka )       = 0._wp
325            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
326            u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)             
327         END DO 
328         !!
329         !! Matrix inversion
330         !! ----------------------------------------------------------
331         DO ji = 2, jpi         
332            zcff                 =   1._wp / z_elem_b( ji, 2 )
333            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 ) 
334            u_abl (ji,jj,2,nt_a) =    zcff * u_abl(ji,jj,2,nt_a)
335         END DO             
336
337         DO jk = 3, jpka           
338            DO ji = 2, jpi             
339               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
340               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
341               u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk  ,nt_a)   &
342               &          - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
343            END DO
344         END DO
345           
346         DO jk = jpkam1,2,-1           
347            DO ji = 2, jpi
348               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)
349            END DO
350         END DO
351                                         
352      !-------------
353      END DO             ! end outer loop
354      !-------------   
355
356      !
357      !  Vertical diffusion for v_abl     
358      !-------------
359      DO jj = 2, jpj   ! outer loop
360      !-------------   
361         !
362         DO jk = 3, jpkam1
363            DO ji = 1, jpi   
364               z_elem_a( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
365               z_elem_c( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal             
366               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal
367            END DO
368         END DO
369
370         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)         
371            !++ Surface boundary condition
372            z_elem_a( ji,     2    ) = 0._wp
373            z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )       
374            !
375         zztmp1 = pcd_du(ji, jj)
376         zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) 
377#if defined key_si3             
378            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
379         zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) )
380         zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
381#endif         
382            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 
383            v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
384            !++ Top Neumann B.C.           
385            !z_elem_a( ji,     jpka ) = -rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
386            !z_elem_c( ji,     jpka ) = 0._wp
387            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
388            !++ Top Dirichlet B.C.
389            z_elem_a( ji,     jpka )       = 0._wp
390            z_elem_c( ji,     jpka )       = 0._wp
391            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
392            v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)             
393         END DO 
394         !!
395         !! Matrix inversion
396         !! ----------------------------------------------------------
397         DO ji = 1, jpi             
398            zcff                 =  1._wp / z_elem_b( ji, 2 )
399            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       ) 
400            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a ) 
401         END DO             
402
403         DO jk = 3, jpka           
404            DO ji = 1, jpi                   
405               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
406               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
407               v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk  ,nt_a)   &
408               &          - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
409            END DO
410         END DO
411           
412         DO jk = jpkam1,2,-1           
413            DO ji = 1, jpi   
414               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)
415            END DO
416         END DO
417         !                                         
418      !-------------
419      END DO             ! end outer loop
420      !-------------   
421   
422      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
423      !                            !  5 *** Apply nudging on the dynamics and the tracers
424      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
425      z_cft(:,:,:) = 0._wp     
426     
427      IF( nn_dyn_restore > 0  ) THEN
428         !-------------
429         DO jk = 2, jpka    ! outer loop
430         !-------------       
431            DO_2D_01_01
432               zcff1 = pblh( ji, jj )
433               zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )                       
434               zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
435               zmsk  = msk_abl(ji,jj)
436               zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   &
437                  &  + jp_alp1_dyn * zsig    + jp_alp0_dyn
438               zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points
439                                                             ! rn_Dt = rDt_abl / nn_fsbc                         
440               zcff  = zcff * rest_eq(ji,jj)
441               z_cft( ji, jj, jk ) = zcff
442               u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   &
443                  &                               + zcff   * pu_dta( ji, jj, jk       )                     
444               v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   &
445                  &                               + zcff   * pv_dta( ji, jj, jk       )
446            END_2D
447         !-------------
448         END DO             ! end outer loop
449         !-------------               
450      END IF
451
452      !-------------
453      DO jk = 2, jpka    ! outer loop
454      !-------------       
455         DO_2D_11_11
456            zcff1 = pblh( ji, jj )
457            zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
458            zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
459            zmsk  = msk_abl(ji,jj)
460            zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   &
461               &  + jp_alp1_tra * zsig    + jp_alp0_tra
462            zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points
463                                                          ! rn_Dt = rDt_abl / nn_fsbc                         
464            !z_cft( ji, jj, jk ) = zcff
465            tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   &
466               &                                       + zcff   * pt_dta( ji, jj, jk )
467           
468            tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa )   &
469               &                                       + zcff   * pq_dta( ji, jj, jk )
470           
471         END_2D
472      !-------------
473      END DO             ! end outer loop
474      !-------------     
475      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
476      !                            !  6 *** MPI exchanges 
477      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
478      !
479      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a      ), 'T', -1. )
480      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...
481      !
482      ! first ABL level
483      IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl",   u_abl(:,:,2,nt_a      ) )
484      IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl",   v_abl(:,:,2,nt_a      ) )
485      IF ( iom_use("tz1_abl") ) CALL iom_put ( "tz1_abl",  tq_abl(:,:,2,nt_a,jp_ta) )
486      IF ( iom_use("qz1_abl") ) CALL iom_put ( "qz1_abl",  tq_abl(:,:,2,nt_a,jp_qa) )
487      IF ( iom_use("uz1_dta") ) CALL iom_put ( "uz1_dta",  pu_dta(:,:,2           ) )
488      IF ( iom_use("vz1_dta") ) CALL iom_put ( "vz1_dta",  pv_dta(:,:,2           ) )
489      IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta",  pt_dta(:,:,2           ) )
490      IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta",  pq_dta(:,:,2           ) )
491      ! all ABL levels
492      IF ( iom_use("u_abl"  ) ) CALL iom_put ( "u_abl"  ,   u_abl(:,:,2:jpka,nt_a      ) )
493      IF ( iom_use("v_abl"  ) ) CALL iom_put ( "v_abl"  ,   v_abl(:,:,2:jpka,nt_a      ) )
494      IF ( iom_use("t_abl"  ) ) CALL iom_put ( "t_abl"  ,  tq_abl(:,:,2:jpka,nt_a,jp_ta) )
495      IF ( iom_use("q_abl"  ) ) CALL iom_put ( "q_abl"  ,  tq_abl(:,:,2:jpka,nt_a,jp_qa) )
496      IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a      ) )
497      IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka           ) )
498      IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka           ) )
499      IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxl_abl(:,:,2:jpka           ) )
500      IF ( iom_use("pblh"   ) ) CALL iom_put (  "pblh"  ,    pblh(:,:                  ) )
501      ! debug (to be removed)
502      IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta",  pu_dta(:,:,2:jpka) )
503      IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta",  pv_dta(:,:,2:jpka) )
504      IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta",  pt_dta(:,:,2:jpka) )
505      IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta",  pq_dta(:,:,2:jpka) )
506      IF ( iom_use("coeft") ) CALL iom_put ( "coeft",   z_cft(:,:,2:jpka) )
507      IF( ln_geos_winds ) THEN
508         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2           ) )
509         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2           ) )
510      END IF
511      IF( ln_hpgls_frc ) THEN
512         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo",  pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)   )
513         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)   )
514      END IF
515      !
516      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
517      !                            !  7 *** Finalize flux computation
518      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
519
520      DO_2D_11_11
521         ztemp             = tq_abl  ( ji, jj, 2, nt_a, jp_ta ) 
522         zhumi             = tq_abl  ( ji, jj, 2, nt_a, jp_qa ) 
523         !zcff              = pslp_dta( ji, jj ) /   &              !<-- At this point ztemp and zhumi should not be zero ...
524         !   &                        (  R_dry*ztemp * ( 1._wp + rctv0*zhumi )  )
525         zcff              = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
526         psen ( ji, jj )   =      cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )
527         pevp ( ji, jj )   = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) )
528         rhoa( ji, jj )   = zcff             
529      END_2D
530     
531      DO_2D_01_01
532         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji  ,jj) + pssu(ji-1,jj) ) 
533         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj  ) + pssv(ji,jj-1) ) 
534      END_2D
535      !
536      CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
537      !
538      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
539      DO_2D_11_11
540         zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   &
541            &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj)  )  ! * msk_abl(ji,jj)
542         zztmp         = rhoa(ji,jj) * pcd_du(ji,jj)
543         
544         pwndm (ji,jj) =         zcff
545         ptaum (ji,jj) = zztmp * zcff
546         zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
547         zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
548      END_2D
549      ! ... utau, vtau at U- and V_points, resp.
550      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
551      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
552      DO_2D_00_00
553         zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
554         zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
555         ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) )
556         zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
557         zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
558         ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) )
559      END_2D
560      !
561      CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1., ptauj(:,:), 'V', -1. )
562
563      CALL iom_put( "taum_oce", ptaum )
564
565      IF(sn_cfctl%l_prtctl) THEN
566         CALL prt_ctl( tab2d_1=pwndm  , clinfo1=' abl_stp: wndm   : ' )
567         CALL prt_ctl( tab2d_1=ptaui  , clinfo1=' abl_stp: utau   : ' )
568         CALL prt_ctl( tab2d_2=ptauj  , clinfo2=          'vtau   : ' )
569      ENDIF
570
571#if defined key_si3
572         ! ------------------------------------------------------------ !
573         !    Wind stress relative to the moving ice ( U10m - U_ice )   !
574         ! ------------------------------------------------------------ !
575         DO_2D_00_00
576           
577            zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) )
578            zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
579   
580            ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             &
581               &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          &
582               &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) )
583            ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             &
584               &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          &
585               &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) )
586         END_2D
587         CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. )
588         !
589         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   &
590            &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' )
591#endif
592      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
593      !                            !  8 *** Swap time indices for the next timestep
594      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
595      nt_n = 1 + MOD( kt  , 2)
596      nt_a = 1 + MOD( kt+1, 2)
597      !   
598!---------------------------------------------------------------------------------------------------
599   END SUBROUTINE abl_stp
600!===================================================================================================
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618!===================================================================================================
619   SUBROUTINE abl_zdf_tke( )
620!---------------------------------------------------------------------------------------------------
621
622      !!----------------------------------------------------------------------
623      !!                   ***  ROUTINE abl_zdf_tke  ***
624      !!
625      !! ** Purpose :   Time-step Turbulente Kinetic Energy (TKE) equation
626      !!
627      !! ** Method  : - source term due to shear
628      !!              - source term due to stratification
629      !!              - resolution of the TKE equation by inverting
630      !!                a tridiagonal linear system
631      !!
632      !! ** Action  : - en : now turbulent kinetic energy)
633      !!              - avmu, avmv : production of TKE by shear at u and v-points
634      !!                (= Kz dz[Ub] * dz[Un] )
635      !! ---------------------------------------------------------------------
636      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn 
637      INTEGER, DIMENSION(1:jpi          )     ::   ikbl
638      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
639      REAL(wp)                                ::   zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2
640      REAL(wp)                                ::   zdU2,zdV2     
641      REAL(wp)                                ::   zwndi,zwndj
642      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2
643      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2
644      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF       
645      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a
646      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b
647      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_c
648      LOGICAL                                 ::   ln_Patankar    = .FALSE.
649      LOGICAL                                 ::   ln_dumpvar     = .FALSE.
650      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl 
651      !
652      tind  = nt_n
653      ziRic = 1._wp / rn_Ric
654      ! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
655      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
656      !                            !  Advance TKE equation to time n+1
657      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
658      !-------------
659      DO jj = 1, jpj    ! outer loop
660      !-------------
661         !
662         ! Compute vertical shear         
663         DO jk = 2, jpkam1
664            DO ji = 1,jpi 
665               zcff        = 1.0_wp / e3w_abl( jk )**2               
666               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2                           
667               zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
668               zsh2(ji,jk) = zdU+zdV
669            END DO
670         END DO   
671         !
672         ! Compute brunt-vaisala frequency
673         DO jk = 2, jpkam1
674            DO ji = 1,jpi 
675               zcff  = grav * itvref / e3w_abl( jk ) 
676               zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta)
677               zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        &
678                  &   - tq_abl( ji, jj, jk  , tind, jp_ta) * tq_abl( ji, jj, jk  , tind, jp_qa)
679               zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi)
680            END DO
681         END DO 
682         !
683         ! Terms for the tridiagonal problem
684         DO jk = 2, jpkam1
685            DO ji = 1,jpi 
686               zshear       =                 zsh2( ji,     jk )   ! zsh2 is already multiplied by Avm_abl at this point
687               zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation         
688               zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
689               
690               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
691               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         
692               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE
693                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
694                     &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal       
695                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )             ! right-hand-side
696               ELSE
697                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
698                     &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal   
699                     &                     - e3w_abl(jk) * rDt_abl * zbuoy   
700                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )             ! right-hand-side                     
701               END IF
702            END DO
703         END DO 
704                           
705         DO ji = 1,jpi    ! vector opt.           
706            zesrf   =  MAX( 4.63_wp * ustar2(ji,jj), tke_min )           
707            zetop   = tke_min             
708            z_elem_a ( ji,     1    ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp                                                   
709            z_elem_a ( ji,     jpka ) = 0._wp; z_elem_c ( ji,     jpka ) = 0._wp; z_elem_b ( ji,     jpka ) = 1._wp           
710            tke_abl( ji, jj,    1, nt_a ) = zesrf
711            tke_abl( ji, jj, jpka, nt_a ) = zetop 
712            zbn2(ji,jj,   1) = zbn2( ji,jj,   2) 
713            zsh2(ji,      1) = zsh2( ji,      2)                                     
714            zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 
715            zsh2(ji,   jpka) = zsh2( ji  , jpkam1)
716         END DO       
717         !!
718         !! Matrix inversion
719         !! ----------------------------------------------------------
720         DO ji = 1,jpi
721            zcff                  =  1._wp / z_elem_b( ji, 1 )
722            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       ) 
723            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a ) 
724         END DO             
725
726         DO jk = 2, jpka           
727            DO ji = 1,jpi
728               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
729               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
730               tke_abl(ji,jj,jk,nt_a) =   zcff * ( tke_abl(ji,jj,jk  ,nt_a)   &
731               &          - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
732            END DO
733         END DO
734           
735         DO jk = jpkam1,1,-1           
736            DO ji = 1,jpi
737               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)
738            END DO
739         END DO
740         
741!!FL should not be needed because of Patankar procedure 
742         tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
743
744         !!
745         !! Diagnose PBL height
746         !! ----------------------------------------------------------
747         
748         
749         !                                                       
750         ! arrays zRH, zFC and zCF are available at this point
751         ! and zFC(:, 1 ) = 0.
752         ! diagnose PBL height based on zsh2 and zbn2
753         zFC (  :  ,1) = 0._wp
754         ikbl( 1:jpi ) = 0 
755         
756         DO jk = 2,jpka
757            DO ji = 1, jpi 
758               zcff  = ghw_abl( jk-1 )
759               zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
760               zcff  = ghw_abl( jk   )
761               zcff2 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
762               zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*(                 &
763                               zcff2 * ( zsh2( ji, jk  ) - ziRic * zbn2( ji, jj, jk   ) &
764                           - rn_Cek  * ( fft_abl( ji, jj  ) * fft_abl( ji, jj ) ) ) &
765                             + zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
766                           - rn_Cek  * ( fft_abl( ji, jj  ) * fft_abl( ji, jj ) ) ) &
767                           &                                                 )
768               IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
769            END DO
770         END DO
771         !
772         ! finalize the computation of the PBL height
773         DO ji = 1, jpi
774            jk = ikbl(ji)
775            IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
776               pblh( ji, jj ) =  (  ghw_abl( jk-1 ) * zFC( ji, jk   )       &
777                  &              -  ghw_abl( jk   ) * zFC( ji, jk-1 )       &
778                  &              ) / ( zFC( ji, jk   ) - zFC( ji, jk-1 ) )
779            ELSE IF( jk==2 ) THEN
780               pblh( ji, jj ) = ghw_abl(2   )
781            ELSE
782               pblh( ji, jj ) = ghw_abl(jpka)
783            END IF
784         END DO
785      !-------------     
786      END DO     
787      !-------------
788      !
789      ! Optional : could add pblh smoothing if pblh is noisy horizontally ...
790      IF(ln_smth_pblh) THEN
791         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
792         CALL smooth_pblh( pblh, msk_abl )
793         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)   
794      ENDIF
795      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
796      !                            !  Diagnostic mixing length computation
797      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
798      !
799      SELECT CASE ( nn_amxl )
800      !
801      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom         
802#   define zlup zRH     
803#   define zldw zFC 
804         DO jj = 1, jpj     ! outer loop
805            !
806            DO ji = 1, jpi
807               mxl_abl ( ji, jj,    1 )  = 0._wp
808               mxl_abl ( ji, jj, jpka )  = mxl_min
809               zldw( ji,        1 )  = 0._wp             
810               zlup( ji,     jpka )  = 0._wp
811            END DO
812            !
813            DO jk = 2, jpkam1   
814               DO ji = 1, jpi 
815                  zbuoy             = MAX( zbn2(ji, jj, jk), rsmall ) 
816                  mxl_abl( ji, jj, jk ) = MAX( mxl_min,  &
817                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   ) 
818               END DO
819            END DO   
820            !
821            ! Limit mxl
822            DO jk = jpkam1,1,-1   
823               DO ji = 1, jpi 
824                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) )
825               END DO
826            END DO   
827            !
828            DO jk = 2, jpka
829               DO ji = 1, jpi 
830                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )   
831               END DO
832            END DO 
833            !
834            DO jk = 1, jpka
835               DO ji = 1, jpi
836                  mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )   
837               END DO
838            END DO                         
839            !
840         END DO
841#   undef zlup     
842#   undef zldw 
843         !
844         !
845      CASE ( 1 )             ! length-scale computed as the distance to the PBL height
846         DO jj = 1,jpj      ! outer loop
847            !
848            DO ji = 1, jpi   ! vector opt.
849               zcff = 1._wp / pblh( ji, jj )     ! inverse of hbl             
850               DO jk = 1, jpka             
851                  zsig  = MIN( zcff * ghw_abl( jk ), 1. )   
852                  zcff1 = pblh( ji, jj )                 
853                  mxl_abl( ji, jj, jk ) =  mxl_min                           &
854                     &  +  zsig         * (  amx1*zcff1 + bmx1*mxl_min ) &
855                     &  +  zsig *  zsig * (  amx2*zcff1 + bmx2*mxl_min ) &
856                     &  +  zsig**3      * (  amx3*zcff1 + bmx3*mxl_min ) &
857                     &  +  zsig**4      * (  amx4*zcff1 + bmx4*mxl_min ) &
858                     &  +  zsig**5      * (  amx5*zcff1 + bmx5*mxl_min )
859               END DO
860            END DO 
861            !           
862         END DO           
863         !
864      CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale
865         !
866#   define zlup zRH     
867#   define zldw zFC           
868! zCF is used for matrix inversion
869!             
870       DO jj = 1, jpj      ! outer loop
871         
872         DO ji = 1, jpi           
873            zlup( ji,    1 ) = mxl_min
874            zldw( ji,    1 ) = mxl_min               
875            zlup( ji, jpka ) = mxl_min
876            zldw( ji, jpka ) = mxl_min               
877         END DO           
878         
879         DO jk = 2,jpka-1
880            DO ji = 1, jpi
881               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
882               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)       
883            END DO
884         END DO         
885         !!
886         !! BL89 search for lup
887         !! ----------------------------------------------------------           
888         DO jk=2,jpka-1 
889            !
890            DO ji = 1, jpi
891               zCF(ji,1:jpka) = 0._wp
892               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
893               ln_foundl(ji ) = .false.
894            END DO   
895            !           
896            DO jkup=jk+1,jpka-1
897               DO ji = 1, jpi
898                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
899                     &               ( zbn2(ji,jj,jkup  )*(ghw_abl(jkup  )-ghw_abl(jk)) &
900                     &               + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )
901                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
902                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
903                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)                   
904                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
905                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
906                     zlup(ji,jk) = zcff               
907                     ln_foundl(ji) = .true.
908                  END IF
909               END DO
910            END DO
911            !
912         END DO   
913         !!
914         !! BL89 search for ldwn
915         !! ----------------------------------------------------------         
916         DO jk=2,jpka-1         
917            !
918            DO ji = 1, jpi
919               zCF(ji,1:jpka) = 0._wp
920               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
921               ln_foundl(ji ) = .false.
922            END DO 
923            !   
924            DO jkdwn=jk-1,1,-1
925               DO ji = 1, jpi             
926                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
927                     &               * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
928                                       + zbn2(ji,jj,jkdwn  )*(ghw_abl(jk)-ghw_abl(jkdwn  )) )   
929                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
930                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
931                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )             
932                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
933                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
934                     zldw(ji,jk) = zcff 
935                     ln_foundl(ji) = .true.               
936                  END IF                                   
937               END DO
938            END DO
939            !     
940         END DO
941
942         DO jk = 1, jpka
943            DO ji = 1, jpi 
944               mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min ) 
945            END DO
946         END DO 
947
948      END DO
949#   undef zlup     
950#   undef zldw           
951         !
952      END SELECT     
953      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
954      !                            !  Finalize the computation of turbulent visc./diff.
955      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
956     
957      !-------------
958      DO jj = 1, jpj     ! outer loop
959      !-------------
960         DO jk = 1, jpka   
961            DO ji = 1, jpi  ! vector opt.
962               zcff              = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk )  &
963                  &                                       * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) ) 
964               zcff2             =  1. / ( 1. + zcff )   !<-- phi_z(z)
965               zcff              = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
966!!FL: MAX function probably useless because of the definition of mxl_min             
967               Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   )
968               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )                             
969            END DO
970         END DO
971      !-------------         
972      END DO     
973      !-------------
974
975!---------------------------------------------------------------------------------------------------
976   END SUBROUTINE abl_zdf_tke
977!===================================================================================================
978
979
980!===================================================================================================
981   SUBROUTINE smooth_pblh( pvar2d, msk )
982!---------------------------------------------------------------------------------------------------
983
984      !!----------------------------------------------------------------------
985      !!                   ***  ROUTINE smooth_pblh  ***
986      !!
987      !! ** Purpose :   2D Hanning filter on atmospheric PBL height
988      !!
989      !! ---------------------------------------------------------------------
990     REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk   
991     REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
992      INTEGER                                     :: ji,jj
993     REAL(wp)                                    :: smth_a, smth_b
994     REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY
995     REAL(wp)                                    :: zumsk,zvmsk
996      !!
997      !!=========================================================
998      !!
999      !! Hanning filter
1000      smth_a = 1._wp / 8._wp
1001      smth_b = 1._wp / 4._wp
1002      !
1003      DO_2D_11_10
1004         zumsk = msk(ji,jj) * msk(ji+1,jj)
1005         zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk
1006      END_2D
1007     
1008     DO_2D_10_11
1009         zvmsk = msk(ji,jj) * msk(ji,jj+1)
1010         zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk
1011     END_2D
1012     
1013     DO_2D_10_00
1014         zFY ( ji, jj  ) =   zdY ( ji, jj   )                        &
1015            & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   &
1016            &            -  (zdX ( ji, jj   ) - zdX( ji-1, jj   ))  )
1017     END_2D
1018
1019      DO_2D_00_10
1020         zFX( ji, jj  ) =    zdX( ji, jj   )                         &
1021           &    + smth_a*(  (zdY( ji+1, jj ) - zdY( ji+1, jj-1))     &
1022           &             -  (zdY( ji  , jj ) - zdY( ji  , jj-1)) )
1023      END_2D
1024
1025      DO_2D_00_00
1026         pvar2d( ji  ,jj ) = pvar2d( ji  ,jj )              &
1027  &         + msk(ji,jj) * smth_b * (                       &
1028  &                  zFX( ji, jj ) - zFX( ji-1, jj )        &
1029  &                 +zFY( ji, jj ) - zFY( ji, jj-1 )  )
1030      END_2D
1031     !!
1032!---------------------------------------------------------------------------------------------------
1033   END SUBROUTINE smooth_pblh
1034!===================================================================================================
1035
1036!!======================================================================
1037END MODULE ablmod
Note: See TracBrowser for help on using the repository browser.