New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
ablmod.F90 in NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL – NEMO

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

Last change on this file since 11937 was 11937, checked in by gsamson, 4 years ago

dev_r11265_ABL: put back lbc_lnk for tracers to ensure mpp reproductibility (#2131)

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