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.
zdfdrg.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfdrg.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • Property svn:keywords set to Id
File size: 25.0 KB
Line 
1MODULE zdfdrg
2   !!======================================================================
3   !!                       ***  MODULE  zdfdrg  ***
4   !! Ocean physics: top and/or Bottom friction
5   !!======================================================================
6   !! History :  OPA  ! 1997-06  (G. Madec, A.-M. Treguier)  Original code
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
8   !!            3.2  ! 2009-09  (A.C.Coward)  Correction to include barotropic contribution
9   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
10   !!            3.4  ! 2011-11  (H. Liu) implementation of semi-implicit bottom friction option
11   !!                 ! 2012-06  (H. Liu) implementation of Log Layer bottom friction option
12   !!            4.0  ! 2017-05  (G. Madec) zdfbfr becomes zdfdrg + variable names change
13   !!                                     + drag defined at t-point + new user interface + top drag (ocean cavities)
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   zdf_drg       : update bottom friction coefficient (non-linear bottom friction only)
18   !!   zdf_drg_exp   : compute the top & bottom friction in explicit case
19   !!   zdf_drg_init  : read in namdrg namelist and control the bottom friction parameters.
20   !!       drg_init  :
21   !!----------------------------------------------------------------------
22   USE oce            ! ocean dynamics and tracers variables
23   USE phycst  , ONLY : vkarmn
24   USE dom_oce        ! ocean space and time domain variables
25   USE zdf_oce        ! ocean vertical physics variables
26   USE trd_oce        ! trends: ocean variables
27   USE trddyn         ! trend manager: dynamics
28   !
29   USE in_out_manager ! I/O manager
30   USE iom            ! I/O module
31   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
32   USE lib_mpp        ! distributed memory computing
33   USE prtctl         ! Print control
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   zdf_drg         ! called by zdf_phy
39   PUBLIC   zdf_drg_exp     ! called by dyn_zdf
40   PUBLIC   zdf_drg_init    ! called by zdf_phy_init
41
42   !                                 !!* Namelist namdrg: nature of drag coefficient namelist *
43   LOGICAL          ::   ln_OFF       ! free-slip       : Cd = 0
44   LOGICAL          ::   ln_lin       !     linear  drag: Cd = Cd0_lin
45   LOGICAL          ::   ln_non_lin   ! non-linear  drag: Cd = Cd0_nl |U|
46   LOGICAL          ::   ln_loglayer  ! logarithmic drag: Cd = vkarmn/log(z/z0)
47   LOGICAL , PUBLIC ::   ln_drgimp    ! implicit top/bottom friction flag
48
49   !                                 !!* Namelist namdrg_top & _bot: TOP or BOTTOM coefficient namelist *
50   REAL(wp)         ::   rn_Cd0       !: drag coefficient                                           [ - ]
51   REAL(wp)         ::   rn_Uc0       !: characteristic velocity (linear case: tau=rho*Cd0*Uc0*u)   [m/s]
52   REAL(wp)         ::   rn_Cdmax     !: drag value maximum (ln_loglayer=T)                         [ - ]
53   REAL(wp)         ::   rn_z0        !: roughness          (ln_loglayer=T)                         [ m ]
54   REAL(wp)         ::   rn_ke0       !: background kinetic energy (non-linear case)                [m2/s2]
55   LOGICAL          ::   ln_boost     !: =T regional boost of Cd0 ; =F Cd0 horizontally uniform
56   REAL(wp)         ::   rn_boost     !: local boost factor                                         [ - ]
57
58   REAL(wp), PUBLIC ::   r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top   ! set from namdrg_top namelist values
59   REAL(wp), PUBLIC ::   r_Cdmin_bot, r_Cdmax_bot, r_z0_bot, r_ke0_bot   !  -    -  namdrg_bot    -       -
60
61   INTEGER ::              ndrg       ! choice of the type of drag coefficient
62   !                                  ! associated indices:
63   INTEGER, PARAMETER ::   np_OFF      = 0   ! free-slip: drag set to zero
64   INTEGER, PARAMETER ::   np_lin      = 1   !     linear drag: Cd = Cd0_lin
65   INTEGER, PARAMETER ::   np_non_lin  = 2   ! non-linear drag: Cd = Cd0_nl |U|
66   INTEGER, PARAMETER ::   np_loglayer = 3   ! non linear drag (logarithmic formulation): Cd = vkarmn/log(z/z0)
67
68   LOGICAL , PUBLIC ::   l_zdfdrg           !: flag to update at each time step the top/bottom Cd
69   LOGICAL          ::   l_log_not_linssh   !: flag to update at each time step the position ot the velocity point
70   !
71   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   rCd0_top, rCd0_bot   !: precomputed top/bottom drag coeff. at t-point (>0)
72   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   rCdU_top, rCdU_bot   !: top/bottom drag coeff. at t-point (<0)  [m/s]
73
74   !! * Substitutions
75#  include "vectopt_loop_substitute.h90"
76   !!----------------------------------------------------------------------
77   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
78   !! $Id$
79   !! Software governed by the CeCILL license (see ./LICENSE)
80   !!----------------------------------------------------------------------
81CONTAINS
82
83   SUBROUTINE zdf_drg( kt, k_mk, pCdmin, pCdmax, pz0, pke0, pCd0,   &   ! <<== in
84      &                                                     pCdU )      ! ==>> out : bottom drag [m/s]
85      !!----------------------------------------------------------------------
86      !!                   ***  ROUTINE zdf_drg  ***
87      !!
88      !! ** Purpose :   update the top/bottom drag coefficient (non-linear case only)
89      !!
90      !! ** Method  :   In non linear friction case, the drag coeficient is
91      !!              a function of the velocity:
92      !!                          Cd = cd0 * |U+Ut|   
93      !!              where U is the top or bottom velocity and
94      !!                    Ut a tidal velocity (Ut^2 = Tidal kinetic energy
95      !!                       assumed here here to be constant)
96      !!              Depending on the input variable, the top- or bottom drag is compted
97      !!
98      !! ** Action  :   p_Cd   drag coefficient at t-point
99      !!----------------------------------------------------------------------
100      INTEGER                 , INTENT(in   ) ::   kt       ! ocean time-step index
101      !                       !               !!         !==  top or bottom variables  ==!
102      INTEGER , DIMENSION(:,:), INTENT(in   ) ::   k_mk     ! wet level (1st or last)
103      REAL(wp)                , INTENT(in   ) ::   pCdmin   ! min drag value
104      REAL(wp)                , INTENT(in   ) ::   pCdmax   ! max drag value
105      REAL(wp)                , INTENT(in   ) ::   pz0      ! roughness
106      REAL(wp)                , INTENT(in   ) ::   pke0     ! background tidal KE
107      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pCd0     ! masked precomputed part of Cd0
108      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCdU     ! = - Cd*|U|   (t-points) [m/s]
109      !!
110      INTEGER ::   ji, jj   ! dummy loop indices
111      INTEGER ::   imk      ! local integers
112      REAL(wp)::   zzz, zut, zvt, zcd   ! local scalars
113      !!----------------------------------------------------------------------
114      !
115      IF( l_log_not_linssh ) THEN     !==  "log layer"  ==!   compute Cd and -Cd*|U|
116         DO jj = 2, jpjm1
117            DO ji = 2, jpim1
118               imk = k_mk(ji,jj)          ! ocean bottom level at t-points
119               zut = un(ji,jj,imk) + un(ji-1,jj,imk)     ! 2 x velocity at t-point
120               zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk)
121               zzz = 0.5_wp * e3t_n(ji,jj,imk)           ! altitude below/above (top/bottom) the boundary
122               !
123!!JC: possible WAD implementation should modify line below if layers vanish
124               zcd = (  vkarmn / LOG( zzz / pz0 )  )**2
125               zcd = pCd0(ji,jj) * MIN(  MAX( pCdmin , zcd ) , pCdmax  )   ! here pCd0 = mask*boost
126               pCdU(ji,jj) = - zcd * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  )
127            END DO
128         END DO
129      ELSE                                            !==  standard Cd  ==!
130         DO jj = 2, jpjm1
131            DO ji = 2, jpim1
132               imk = k_mk(ji,jj)    ! ocean bottom level at t-points
133               zut = un(ji,jj,imk) + un(ji-1,jj,imk)     ! 2 x velocity at t-point
134               zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk)
135               !                                                           ! here pCd0 = mask*boost * drag
136               pCdU(ji,jj) = - pCd0(ji,jj) * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  )
137            END DO
138         END DO
139      ENDIF
140      !
141      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pCdU, clinfo1=' Cd*U ')
142      !
143   END SUBROUTINE zdf_drg
144
145
146   SUBROUTINE zdf_drg_exp( kt, Kmm, pub, pvb, pua, pva )
147      !!----------------------------------------------------------------------
148      !!                  ***  ROUTINE zdf_drg_exp  ***
149      !!
150      !! ** Purpose :   compute and add the explicit top and bottom frictions.
151      !!
152      !! ** Method  :   in explicit case,
153      !!
154      !!              NB: in implicit case the calculation is performed in dynzdf.F90
155      !!
156      !! ** Action  :   (pua,pva)   momentum trend increased by top & bottom friction trend
157      !!---------------------------------------------------------------------
158      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
159      INTEGER                         , INTENT(in   ) ::   Kmm        ! time level indices
160      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pub, pvb   ! the two components of the before velocity
161      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! the two components of the velocity tendency
162      !!
163      INTEGER  ::   ji, jj       ! dummy loop indexes
164      INTEGER  ::   ikbu, ikbv   ! local integers
165      REAL(wp) ::   zm1_2dt      ! local scalar
166      REAL(wp) ::   zCdu, zCdv   !   -      -
167      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv
168      !!---------------------------------------------------------------------
169      !
170!!gm bug : time step is only rdt (not 2 rdt if euler start !)
171      zm1_2dt = - 1._wp / ( 2._wp * rdt )
172
173      IF( l_trddyn ) THEN      ! trends: store the input trends
174         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
175         ztrdu(:,:,:) = pua(:,:,:)
176         ztrdv(:,:,:) = pva(:,:,:)
177      ENDIF
178
179      DO jj = 2, jpjm1
180         DO ji = 2, jpim1
181            ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
182            ikbv = mbkv(ji,jj)
183            !
184            ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
185            zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu)
186            zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv)
187            !
188            pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
189            pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv)
190         END DO
191      END DO
192      !
193      IF( ln_isfcav ) THEN        ! ocean cavities
194         DO jj = 2, jpjm1
195            DO ji = 2, jpim1
196               ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
197               ikbv = mikv(ji,jj)
198               !
199               ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
200               zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu)    ! NB: Cdtop masked
201               zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv)
202               !
203               pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
204               pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv)
205           END DO
206         END DO
207      ENDIF
208      !
209      IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics
210         ztrdu(:,:,:) = pua(:,:,:) - ztrdu(:,:,:)
211         ztrdv(:,:,:) = pva(:,:,:) - ztrdv(:,:,:)
212         CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt, Kmm )
213         DEALLOCATE( ztrdu, ztrdv )
214      ENDIF
215      !                                          ! print mean trends (used for debugging)
216      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pua, clinfo1=' bfr  - Ua: ', mask1=umask,               &
217         &                       tab3d_2=pva, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
218      !
219   END SUBROUTINE zdf_drg_exp
220
221
222   SUBROUTINE zdf_drg_init
223      !!----------------------------------------------------------------------
224      !!                  ***  ROUTINE zdf_brg_init  ***
225      !!
226      !! ** Purpose :   Initialization of the bottom friction
227      !!
228      !! ** Method  :   Read the namdrg namelist and check their consistency
229      !!                called at the first timestep (nit000)
230      !!----------------------------------------------------------------------
231      INTEGER   ::   ji, jj      ! dummy loop indexes
232      INTEGER   ::   ios, ioptio   ! local integers
233      !!
234      NAMELIST/namdrg/ ln_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp
235      !!----------------------------------------------------------------------
236      !
237      !                     !==  drag nature  ==!
238      !
239      REWIND( numnam_ref )                   ! Namelist namdrg in reference namelist
240      READ  ( numnam_ref, namdrg, IOSTAT = ios, ERR = 901)
241901   IF( ios /= 0 )   CALL ctl_nam( ios , 'namdrg in reference namelist', lwp )
242      REWIND( numnam_cfg )                   ! Namelist namdrg in configuration namelist
243      READ  ( numnam_cfg, namdrg, IOSTAT = ios, ERR = 902 )
244902   IF( ios >  0 )   CALL ctl_nam( ios , 'namdrg in configuration namelist', lwp )
245      IF(lwm) WRITE ( numond, namdrg )
246      !
247      IF(lwp) THEN
248         WRITE(numout,*)
249         WRITE(numout,*) 'zdf_drg_init : top and/or bottom drag setting'
250         WRITE(numout,*) '~~~~~~~~~~~~'
251         WRITE(numout,*) '   Namelist namdrg : top/bottom friction choices'
252         WRITE(numout,*) '      free-slip       : Cd = 0                  ln_OFF      = ', ln_OFF 
253         WRITE(numout,*) '      linear  drag    : Cd = Cd0                ln_lin      = ', ln_lin
254         WRITE(numout,*) '      non-linear  drag: Cd = Cd0_nl |U|         ln_non_lin  = ', ln_non_lin
255         WRITE(numout,*) '      logarithmic drag: Cd = vkarmn/log(z/z0)   ln_loglayer = ', ln_loglayer
256         WRITE(numout,*) '      implicit friction                         ln_drgimp   = ', ln_drgimp
257      ENDIF
258      !
259      ioptio = 0                       ! set ndrg and control check
260      IF( ln_OFF      ) THEN   ;   ndrg = np_OFF        ;   ioptio = ioptio + 1   ;   ENDIF
261      IF( ln_lin      ) THEN   ;   ndrg = np_lin        ;   ioptio = ioptio + 1   ;   ENDIF
262      IF( ln_non_lin  ) THEN   ;   ndrg = np_non_lin    ;   ioptio = ioptio + 1   ;   ENDIF
263      IF( ln_loglayer ) THEN   ;   ndrg = np_loglayer   ;   ioptio = ioptio + 1   ;   ENDIF
264      !
265      IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' )
266      !
267      !
268      !                     !==  BOTTOM drag setting  ==!   (applied at seafloor)
269      !
270      ALLOCATE( rCd0_bot(jpi,jpj), rCdU_bot(jpi,jpj) )
271      CALL drg_init( 'BOTTOM'   , mbkt       ,                                         &   ! <== in
272         &           r_Cdmin_bot, r_Cdmax_bot, r_z0_bot, r_ke0_bot, rCd0_bot, rCdU_bot )   ! ==> out
273
274      !
275      !                     !==  TOP drag setting  ==!   (applied at the top of ocean cavities)
276      !
277      IF( ln_isfcav ) THEN              ! Ocean cavities: top friction setting
278         ALLOCATE( rCd0_top(jpi,jpj), rCdU_top(jpi,jpj) )
279         CALL drg_init( 'TOP   '   , mikt       ,                                         &   ! <== in
280            &           r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top, rCd0_top, rCdU_top )   ! ==> out
281      ENDIF
282      !
283   END SUBROUTINE zdf_drg_init
284
285
286   SUBROUTINE drg_init( cd_topbot, k_mk,  &
287      &                 pCdmin, pCdmax, pz0, pke0, pCd0, pCdU ) 
288      !!----------------------------------------------------------------------
289      !!                  ***  ROUTINE drg_init  ***
290      !!
291      !! ** Purpose :   Initialization of the top/bottom friction CdO and Cd
292      !!              from namelist parameters
293      !!----------------------------------------------------------------------
294      CHARACTER(len=6)        , INTENT(in   ) ::   cd_topbot       ! top/ bot indicator
295      INTEGER , DIMENSION(:,:), INTENT(in   ) ::   k_mk            ! 1st/last  wet level
296      REAL(wp)                , INTENT(  out) ::   pCdmin, pCdmax  ! min and max drag coef. [-]
297      REAL(wp)                , INTENT(  out) ::   pz0             ! roughness              [m]
298      REAL(wp)                , INTENT(  out) ::   pke0            ! background KE          [m2/s2]
299      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCd0            ! masked precomputed part of the non-linear drag coefficient
300      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCdU            ! minus linear drag*|U| at t-points  [m/s]
301      !!
302      CHARACTER(len=40) ::   cl_namdrg, cl_file, cl_varname, cl_namref, cl_namcfg  ! local names
303      INTEGER ::   ji, jj              ! dummy loop indexes
304      LOGICAL ::   ll_top, ll_bot      ! local logical
305      INTEGER ::   ios, inum, imk      ! local integers
306      REAL(wp)::   zmsk, zzz, zcd      ! local scalars
307      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk_boost   ! 2D workspace
308      !!
309      NAMELIST/namdrg_top/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost
310      NAMELIST/namdrg_bot/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost
311      !!----------------------------------------------------------------------
312      !
313      !                          !==  set TOP / BOTTOM specificities  ==!
314      ll_top = .FALSE.
315      ll_bot = .FALSE.
316      !
317      SELECT CASE (cd_topbot)
318      CASE( 'TOP   ' )
319         ll_top = .TRUE.
320         cl_namdrg  = 'namdrg_top'
321         cl_namref  = 'namdrg_top in reference     namelist'
322         cl_namcfg  = 'namdrg_top in configuration namelist'
323         cl_file    = 'tfr_coef.nc'
324         cl_varname = 'tfr_coef'
325      CASE( 'BOTTOM' )
326         ll_bot = .TRUE.
327         cl_namdrg  = 'namdrg_bot'
328         cl_namref  = 'namdrg_bot  in reference     namelist'
329         cl_namcfg  = 'namdrg_bot  in configuration namelist'
330         cl_file    = 'bfr_coef.nc'
331         cl_varname = 'bfr_coef'
332      CASE DEFAULT
333         CALL ctl_stop( 'drg_init: bad value for cd_topbot ' )
334      END SELECT
335      !
336      !                          !==  read namlist  ==!
337      !
338      REWIND( numnam_ref )                   ! Namelist cl_namdrg in reference namelist
339      IF(ll_top)   READ  ( numnam_ref, namdrg_top, IOSTAT = ios, ERR = 901)
340      IF(ll_bot)   READ  ( numnam_ref, namdrg_bot, IOSTAT = ios, ERR = 901)
341901   IF( ios /= 0 )   CALL ctl_nam( ios , TRIM(cl_namref), lwp )
342      REWIND( numnam_cfg )                   ! Namelist cd_namdrg in configuration namelist
343      IF(ll_top)   READ  ( numnam_cfg, namdrg_top, IOSTAT = ios, ERR = 902 )
344      IF(ll_bot)   READ  ( numnam_cfg, namdrg_bot, IOSTAT = ios, ERR = 902 )
345902   IF( ios >  0 )   CALL ctl_nam( ios , TRIM(cl_namcfg), lwp )
346      IF(lwm .AND. ll_top)   WRITE ( numond, namdrg_top )
347      IF(lwm .AND. ll_bot)   WRITE ( numond, namdrg_bot )
348      !
349      IF(lwp) THEN
350         WRITE(numout,*)
351         WRITE(numout,*) '   Namelist ',TRIM(cl_namdrg),' : set ',TRIM(cd_topbot),' friction parameters'
352         WRITE(numout,*) '      drag coefficient                        rn_Cd0   = ', rn_Cd0
353         WRITE(numout,*) '      characteristic velocity (linear case)   rn_Uc0   = ', rn_Uc0, ' m/s'
354         WRITE(numout,*) '      non-linear drag maximum                 rn_Cdmax = ', rn_Cdmax
355         WRITE(numout,*) '      background kinetic energy  (n-l case)   rn_ke0   = ', rn_ke0
356         WRITE(numout,*) '      bottom roughness           (n-l case)   rn_z0    = ', rn_z0
357         WRITE(numout,*) '      set a regional boost of Cd0             ln_boost = ', ln_boost
358         WRITE(numout,*) '         associated boost factor              rn_boost = ', rn_boost
359      ENDIF
360      !
361      !                          !==  return some namelist parametres  ==!   (used in non_lin and loglayer cases)
362      pCdmin = rn_Cd0
363      pCdmax = rn_Cdmax
364      pz0    = rn_z0
365      pke0   = rn_ke0
366      !
367      !                          !==  mask * boost factor  ==!
368      !
369      IF( ln_boost ) THEN           !* regional boost:   boost factor = 1 + regional boost
370         IF(lwp) WRITE(numout,*)
371         IF(lwp) WRITE(numout,*) '   ==>>>   use a regional boost read in ', TRIM(cl_file), ' file'
372         IF(lwp) WRITE(numout,*) '           using enhancement factor of ', rn_boost
373         ! cl_varname is a coefficient in [0,1] giving where to apply the regional boost
374         CALL iom_open ( TRIM(cl_file), inum )
375         CALL iom_get  ( inum, jpdom_data, TRIM(cl_varname), zmsk_boost, 1 )
376         CALL iom_close( inum)
377         zmsk_boost(:,:) = 1._wp + rn_boost * zmsk_boost(:,:)
378         !
379      ELSE                          !* no boost:   boost factor = 1
380         zmsk_boost(:,:) = 1._wp
381      ENDIF
382      !                             !* mask outside ocean cavities area (top) or land area (bot)
383      IF(ll_top)   zmsk_boost(:,:) = zmsk_boost(:,:) * ssmask(:,:) * (1. - tmask(:,:,1) )  ! none zero in ocean cavities only
384      IF(ll_bot)   zmsk_boost(:,:) = zmsk_boost(:,:) * ssmask(:,:)                         ! x seafloor mask
385      !
386      !
387      SELECT CASE( ndrg )
388      !
389      CASE( np_OFF  )            !==  No top/bottom friction  ==!   (pCdU = 0)
390         IF(lwp) WRITE(numout,*)
391         IF(lwp) WRITE(numout,*) '   ==>>>   ',TRIM(cd_topbot),' free-slip, friction set to zero'
392         !
393         l_zdfdrg = .FALSE.         ! no time variation of the drag: set it one for all
394         !
395         pCdU(:,:) = 0._wp
396         pCd0(:,:) = 0._wp
397         !
398      CASE( np_lin )             !==  linear friction  ==!   (pCdU = Cd0 * Uc0)
399         IF(lwp) WRITE(numout,*)
400         IF(lwp) WRITE(numout,*) '   ==>>>   linear ',TRIM(cd_topbot),' friction (constant coef = Cd0*Uc0 = ', rn_Cd0*rn_Uc0, ')'
401         !
402         l_zdfdrg = .FALSE.         ! no time variation of the Cd*|U| : set it one for all
403         !                     
404         pCd0(:,:) = rn_Cd0 * zmsk_boost(:,:)  !* constant in time drag coefficient (= mask (and boost) Cd0)
405         pCdU(:,:) = - pCd0(:,:) * rn_Uc0      !  using a constant velocity
406         !
407      CASE( np_non_lin )         !== non-linear friction  ==!   (pCd0 = Cd0 )
408         IF(lwp) WRITE(numout,*)
409         IF(lwp) WRITE(numout,*) '   ==>>>   quadratic ',TRIM(cd_topbot),' friction (propotional to module of the velocity)'
410         IF(lwp) WRITE(numout,*) '   with    a drag coefficient Cd0 = ', rn_Cd0, ', and'
411         IF(lwp) WRITE(numout,*) '           a background velocity module of (rn_ke0)^1/2 = ', SQRT(rn_ke0), 'm/s)'
412         !
413         l_zdfdrg = .TRUE.          !* Cd*|U| updated at each time-step (it depends on ocean velocity)
414         !
415         pCd0(:,:) = rn_Cd0 * zmsk_boost(:,:)  !* constant in time proportionality coefficient (= mask (and boost) Cd0)
416         pCdU(:,:) = 0._wp                     ! 
417         !
418      CASE( np_loglayer )       !== logarithmic layer formulation of friction  ==!   (CdU = (vkarman log(z/z0))^2 |U| )
419         IF(lwp) WRITE(numout,*)
420         IF(lwp) WRITE(numout,*) '   ==>>>   quadratic ',TRIM(cd_topbot),' drag (propotional to module of the velocity)'
421         IF(lwp) WRITE(numout,*) '   with   a logarithmic Cd0 formulation Cd0 = ( vkarman log(z/z0) )^2 ,'
422         IF(lwp) WRITE(numout,*) '          a background velocity module of (rn_ke0)^1/2 = ', SQRT(pke0), 'm/s), '
423         IF(lwp) WRITE(numout,*) '          a logarithmic formulation: a roughness of ', pz0, ' meters,   and '
424         IF(lwp) WRITE(numout,*) '          a proportionality factor bounded by min/max values of ', pCdmin, pCdmax
425         !
426         l_zdfdrg = .TRUE.          !* Cd*|U| updated at each time-step (it depends on ocean velocity)
427         !
428         IF( ln_linssh ) THEN       !* pCd0 = (v log(z/z0))^2   as velocity points have a fixed z position
429            IF(lwp) WRITE(numout,*)
430            IF(lwp) WRITE(numout,*) '   N.B.   linear free surface case, Cd0 computed one for all'
431            !
432            l_log_not_linssh = .FALSE.    !- don't update Cd at each time step
433            !
434            DO jj = 1, jpj                   ! pCd0 = mask (and boosted) logarithmic drag coef.
435               DO ji = 1, jpi
436                  zzz =  0.5_wp * e3t_0(ji,jj,k_mk(ji,jj))
437                  zcd = (  vkarmn / LOG( zzz / rn_z0 )  )**2
438                  pCd0(ji,jj) = zmsk_boost(ji,jj) * MIN(  MAX( rn_Cd0 , zcd ) , rn_Cdmax  )  ! rn_Cd0 < Cd0 < rn_Cdmax
439               END DO
440            END DO
441         ELSE                       !* Cd updated at each time-step ==> pCd0 = mask * boost
442            IF(lwp) WRITE(numout,*)
443            IF(lwp) WRITE(numout,*) '   N.B.   non-linear free surface case, Cd0 updated at each time-step '
444            !
445            l_log_not_linssh = .TRUE.     ! compute the drag coef. at each time-step
446            !
447            pCd0(:,:) = zmsk_boost(:,:)
448         ENDIF
449         pCdU(:,:) = 0._wp          ! initialisation to zero (will be updated at each time step)
450         !
451      CASE DEFAULT
452         CALL ctl_stop( 'drg_init: bad flag value for ndrg ' )
453      END SELECT
454      !
455   END SUBROUTINE drg_init
456
457   !!======================================================================
458END MODULE zdfdrg
Note: See TracBrowser for help on using the repository browser.