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_r11943_MERGE_2019/src/ABL – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ABL/ablmod.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

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