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/CNRS/dev_r14723_tides_under_isf/src/OCE/ZDF – NEMO

source: NEMO/branches/CNRS/dev_r14723_tides_under_isf/src/OCE/ZDF/zdfdrg.F90 @ 15605

Last change on this file since 15605 was 15605, checked in by khutchinson, 3 years ago

changes required to add background tidal velocities in isf cavities (not sette tested yet)

  • Property svn:keywords set to Id
File size: 28.2 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   USE sbc_oce , ONLY : nn_ice 
35! tipaccs 2d top tidal velocity
36   USE fldread, ONLY : FLD_N
37! end tipaccs 2d top tidal velocity
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   zdf_drg         ! called by zdf_phy
43   PUBLIC   zdf_drg_exp     ! called by dyn_zdf
44   PUBLIC   zdf_drg_init    ! called by zdf_phy_init
45
46   !                                 !!* Namelist namdrg: nature of drag coefficient namelist *
47   LOGICAL , PUBLIC ::   ln_drg_OFF   ! free-slip       : Cd = 0
48   LOGICAL          ::   ln_lin       !     linear  drag: Cd = Cd0_lin
49   LOGICAL          ::   ln_non_lin   ! non-linear  drag: Cd = Cd0_nl |U|
50   LOGICAL          ::   ln_loglayer  ! logarithmic drag: Cd = vkarmn/log(z/z0)
51   LOGICAL , PUBLIC ::   ln_drgimp    ! implicit top/bottom friction flag
52   LOGICAL , PUBLIC ::   ln_drgice_imp ! implicit ice-ocean drag
53! tipaccs 2d top tidal velocity
54   LOGICAL , PUBLIC ::   ln_2d_ttv    ! top tidal velocity
55! end tipaccs 2d top tidal velocity
56   !                                 !!* Namelist namdrg_top & _bot: TOP or BOTTOM coefficient namelist *
57   REAL(wp)         ::   rn_Cd0       !: drag coefficient                                           [ - ]
58   REAL(wp)         ::   rn_Uc0       !: characteristic velocity (linear case: tau=rho*Cd0*Uc0*u)   [m/s]
59   REAL(wp)         ::   rn_Cdmax     !: drag value maximum (ln_loglayer=T)                         [ - ]
60   REAL(wp)         ::   rn_z0        !: roughness          (ln_loglayer=T)                         [ m ]
61   REAL(wp)         ::   rn_ke0       !: background kinetic energy (non-linear case)                [m2/s2]
62   LOGICAL          ::   ln_boost     !: =T regional boost of Cd0 ; =F Cd0 horizontally uniform
63   REAL(wp)         ::   rn_boost     !: local boost factor                                         [ - ]
64
65! tipaccs 2d top tidal velocity
66   REAL(wp), PUBLIC ::   r_Cdmin_top, r_Cdmax_top, r_z0_top   ! set from namdrg_top namelist values
67   REAL(wp), PUBLIC ::   r_Cdmin_bot, r_Cdmax_bot, r_z0_bot   !  -    -  namdrg_bot    -       -
68! end tipaccs 2d top tidal velocity
69
70   INTEGER ::              ndrg       ! choice of the type of drag coefficient
71   !                                  ! associated indices:
72   INTEGER, PARAMETER ::   np_OFF      = 0   ! free-slip: drag set to zero
73   INTEGER, PARAMETER ::   np_lin      = 1   !     linear drag: Cd = Cd0_lin
74   INTEGER, PARAMETER ::   np_non_lin  = 2   ! non-linear drag: Cd = Cd0_nl |U|
75   INTEGER, PARAMETER ::   np_loglayer = 3   ! non linear drag (logarithmic formulation): Cd = vkarmn/log(z/z0)
76
77   LOGICAL , PUBLIC ::   l_zdfdrg           !: flag to update at each time step the top/bottom Cd
78   LOGICAL          ::   l_log_not_linssh   !: flag to update at each time step the position ot the velocity point
79   !
80! tipaccs 2d top tidal velocity
81   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   rke0_top, rke0_bot   !: precomputed top/bottom tidal velocity
82! end tipaccs 2d top tidal velocity
83   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   rCd0_top, rCd0_bot   !: precomputed top/bottom drag coeff. at t-point (>0)
84   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   rCdU_top, rCdU_bot   !: top/bottom drag coeff. at t-point (<0)  [m/s]
85
86   !! * Substitutions
87#  include "do_loop_substitute.h90"
88#  include "domzgr_substitute.h90"
89   !!----------------------------------------------------------------------
90   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
91   !! $Id$
92   !! Software governed by the CeCILL license (see ./LICENSE)
93   !!----------------------------------------------------------------------
94CONTAINS
95
96   SUBROUTINE zdf_drg( kt, Kmm, k_mk, pCdmin, pCdmax, pz0, pke0, pCd0,   &   ! <<== in
97      &                                                     pCdU )      ! ==>> out : bottom drag [m/s]
98      !!----------------------------------------------------------------------
99      !!                   ***  ROUTINE zdf_drg  ***
100      !!
101      !! ** Purpose :   update the top/bottom drag coefficient (non-linear case only)
102      !!
103      !! ** Method  :   In non linear friction case, the drag coeficient is
104      !!              a function of the velocity:
105      !!                          Cd = cd0 * |U+Ut|   
106      !!              where U is the top or bottom velocity and
107      !!                    Ut a tidal velocity (Ut^2 = Tidal kinetic energy
108      !!                       assumed here here to be constant)
109      !!              Depending on the input variable, the top- or bottom drag is compted
110      !!
111      !! ** Action  :   p_Cd   drag coefficient at t-point
112      !!----------------------------------------------------------------------
113      INTEGER                 , INTENT(in   ) ::   kt       ! ocean time-step index
114      INTEGER                 , INTENT(in   ) ::   Kmm      ! ocean time level index
115      !                       !               !!         !==  top or bottom variables  ==!
116      INTEGER , DIMENSION(:,:), INTENT(in   ) ::   k_mk     ! wet level (1st or last)
117      REAL(wp)                , INTENT(in   ) ::   pCdmin   ! min drag value
118      REAL(wp)                , INTENT(in   ) ::   pCdmax   ! max drag value
119      REAL(wp)                , INTENT(in   ) ::   pz0      ! roughness
120! tipaccs 2d top tidal velocity
121      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pke0     ! background tidal KE
122! end tipaccs 2d top tidal velocity
123      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pCd0     ! masked precomputed part of Cd0
124      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCdU     ! = - Cd*|U|   (t-points) [m/s]
125      !!
126      INTEGER ::   ji, jj   ! dummy loop indices
127      INTEGER ::   imk      ! local integers
128      REAL(wp)::   zzz, zut, zvt, zcd   ! local scalars
129      !!----------------------------------------------------------------------
130      !
131      IF( l_log_not_linssh ) THEN     !==  "log layer"  ==!   compute Cd and -Cd*|U|
132         DO_2D( 0, 0, 0, 0 )
133            imk = k_mk(ji,jj)          ! ocean bottom level at t-points
134            zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point
135            zvt = vv(ji,jj,imk,Kmm) + vv(ji,jj-1,imk,Kmm)
136            zzz = 0.5_wp * e3t(ji,jj,imk,Kmm)           ! altitude below/above (top/bottom) the boundary
137            !
138!!JC: possible WAD implementation should modify line below if layers vanish
139            zcd = (  vkarmn / LOG( zzz / pz0 )  )**2
140            zcd = pCd0(ji,jj) * MIN(  MAX( pCdmin , zcd ) , pCdmax  )   ! here pCd0 = mask*boost
141            pCdU(ji,jj) = - zcd * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0(ji,jj)  )
142         END_2D
143      ELSE                                            !==  standard Cd  ==!
144         DO_2D( 0, 0, 0, 0 )
145            imk = k_mk(ji,jj)    ! ocean bottom level at t-points
146            zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point
147            zvt = vv(ji,jj,imk,Kmm) + vv(ji,jj-1,imk,Kmm)
148            !                                                           ! here pCd0 = mask*boost * drag
149            pCdU(ji,jj) = - pCd0(ji,jj) * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0(ji,jj)  )
150         END_2D
151      ENDIF
152      !
153      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pCdU, clinfo1=' Cd*U ')
154      !
155   END SUBROUTINE zdf_drg
156
157
158   SUBROUTINE zdf_drg_exp( kt, Kmm, pub, pvb, pua, pva )
159      !!----------------------------------------------------------------------
160      !!                  ***  ROUTINE zdf_drg_exp  ***
161      !!
162      !! ** Purpose :   compute and add the explicit top and bottom frictions.
163      !!
164      !! ** Method  :   in explicit case,
165      !!
166      !!              NB: in implicit case the calculation is performed in dynzdf.F90
167      !!
168      !! ** Action  :   (pua,pva)   momentum trend increased by top & bottom friction trend
169      !!---------------------------------------------------------------------
170      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
171      INTEGER                         , INTENT(in   ) ::   Kmm        ! time level indices
172      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pub, pvb   ! the two components of the before velocity
173      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! the two components of the velocity tendency
174      !!
175      INTEGER  ::   ji, jj       ! dummy loop indexes
176      INTEGER  ::   ikbu, ikbv   ! local integers
177      REAL(wp) ::   zm1_2dt      ! local scalar
178      REAL(wp) ::   zCdu, zCdv   !   -      -
179      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv
180      !!---------------------------------------------------------------------
181      !
182!!gm bug : time step is only rn_Dt (not 2 rn_Dt if euler start !)
183      zm1_2dt = - 1._wp / ( 2._wp * rn_Dt )
184
185      IF( l_trddyn ) THEN      ! trends: store the input trends
186         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
187         ztrdu(:,:,:) = pua(:,:,:)
188         ztrdv(:,:,:) = pva(:,:,:)
189      ENDIF
190
191      DO_2D( 0, 0, 0, 0 )
192         ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
193         ikbv = mbkv(ji,jj)
194         !
195         ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
196         zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u(ji,jj,ikbu,Kmm)
197         zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v(ji,jj,ikbv,Kmm)
198         !
199         pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
200         pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv)
201      END_2D
202      !
203      IF( ln_isfcav ) THEN        ! ocean cavities
204         DO_2D( 0, 0, 0, 0 )
205            ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
206            ikbv = mikv(ji,jj)
207            !
208            ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
209            zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u(ji,jj,ikbu,Kmm)    ! NB: Cdtop masked
210            zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v(ji,jj,ikbv,Kmm)
211            !
212            pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
213            pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv)
214         END_2D
215      ENDIF
216      !
217      IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics
218         ztrdu(:,:,:) = pua(:,:,:) - ztrdu(:,:,:)
219         ztrdv(:,:,:) = pva(:,:,:) - ztrdv(:,:,:)
220         CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt, Kmm )
221         DEALLOCATE( ztrdu, ztrdv )
222      ENDIF
223      !                                          ! print mean trends (used for debugging)
224      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pua, clinfo1=' bfr  - Ua: ', mask1=umask,               &
225         &                                  tab3d_2=pva, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
226      !
227   END SUBROUTINE zdf_drg_exp
228
229
230   SUBROUTINE zdf_drg_init
231      !!----------------------------------------------------------------------
232      !!                  ***  ROUTINE zdf_brg_init  ***
233      !!
234      !! ** Purpose :   Initialization of the bottom friction
235      !!
236      !! ** Method  :   Read the namdrg namelist and check their consistency
237      !!                called at the first timestep (nit000)
238      !!----------------------------------------------------------------------
239      INTEGER   ::   ji, jj      ! dummy loop indexes
240      INTEGER   ::   ios, ioptio   ! local integers
241      !!
242      NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp, ln_drgice_imp
243      !!----------------------------------------------------------------------
244      !
245      !                     !==  drag nature  ==!
246      !
247      READ  ( numnam_ref, namdrg, IOSTAT = ios, ERR = 901)
248901   IF( ios /= 0 )   CALL ctl_nam( ios , 'namdrg in reference namelist' )
249      READ  ( numnam_cfg, namdrg, IOSTAT = ios, ERR = 902 )
250902   IF( ios >  0 )   CALL ctl_nam( ios , 'namdrg in configuration namelist' )
251      IF(lwm) WRITE ( numond, namdrg )
252      !
253      IF ( ln_drgice_imp .AND.   nn_ice /= 2  )   ln_drgice_imp = .FALSE.
254      !
255      IF(lwp) THEN
256         WRITE(numout,*)
257         WRITE(numout,*) 'zdf_drg_init : top and/or bottom drag setting'
258         WRITE(numout,*) '~~~~~~~~~~~~'
259         WRITE(numout,*) '   Namelist namdrg : top/bottom friction choices'
260         WRITE(numout,*) '      free-slip       : Cd = 0                  ln_drg_OFF  = ', ln_drg_OFF 
261         WRITE(numout,*) '      linear  drag    : Cd = Cd0                ln_lin      = ', ln_lin
262         WRITE(numout,*) '      non-linear  drag: Cd = Cd0_nl |U|         ln_non_lin  = ', ln_non_lin
263         WRITE(numout,*) '      logarithmic drag: Cd = vkarmn/log(z/z0)   ln_loglayer = ', ln_loglayer
264         WRITE(numout,*) '      implicit friction                         ln_drgimp   = ', ln_drgimp
265         WRITE(numout,*) '      implicit ice-ocean drag                   ln_drgice_imp  =', ln_drgice_imp
266      ENDIF
267      !
268      ioptio = 0                       ! set ndrg and control check
269      IF( ln_drg_OFF  ) THEN   ;   ndrg = np_OFF        ;   ioptio = ioptio + 1   ;   ENDIF
270      IF( ln_lin      ) THEN   ;   ndrg = np_lin        ;   ioptio = ioptio + 1   ;   ENDIF
271      IF( ln_non_lin  ) THEN   ;   ndrg = np_non_lin    ;   ioptio = ioptio + 1   ;   ENDIF
272      IF( ln_loglayer ) THEN   ;   ndrg = np_loglayer   ;   ioptio = ioptio + 1   ;   ENDIF
273      !
274      IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' )
275      !
276      IF ( ln_drgice_imp.AND.(.NOT.ln_drgimp) ) & 
277         &                CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires ln_drgimp=T' )
278      !
279      !                     !==  BOTTOM drag setting  ==!   (applied at seafloor)
280      !
281      ALLOCATE( rCd0_bot(jpi,jpj), rCdU_bot(jpi,jpj), rke0_bot(jpi,jpj) )
282! tipaccs 2d top tidal velocity
283      CALL drg_init( 'BOTTOM'   , mbkt       ,                                         &   ! <== in
284         &           r_Cdmin_bot, r_Cdmax_bot, r_z0_bot, rke0_bot, rCd0_bot, rCdU_bot )   ! ==> out
285! end tipaccs 2d top tidal velocity
286
287      !
288      !                     !==  TOP drag setting  ==!   (applied at the top of ocean cavities)
289      !
290      IF( ln_isfcav.OR.ln_drgice_imp ) THEN              ! Ocean cavities: top friction setting
291         ALLOCATE( rCdU_top(jpi,jpj) )
292      ENDIF
293      !
294      IF( ln_isfcav ) THEN
295         ALLOCATE( rCd0_top(jpi,jpj), rke0_top(jpi,jpj))
296! tipaccs 2d top tidal velocity
297         CALL drg_init( 'TOP   '   , mikt       ,                                         &   ! <== in
298            &           r_Cdmin_top, r_Cdmax_top, r_z0_top, rke0_top, rCd0_top, rCdU_top )   ! ==> out
299! end tipaccs 2d top tidal velocity
300      ENDIF
301      !
302   END SUBROUTINE zdf_drg_init
303
304
305   SUBROUTINE drg_init( cd_topbot, k_mk,  &
306      &                 pCdmin, pCdmax, pz0, pke0, pCd0, pCdU ) 
307      !!----------------------------------------------------------------------
308      !!                  ***  ROUTINE drg_init  ***
309      !!
310      !! ** Purpose :   Initialization of the top/bottom friction CdO and Cd
311      !!              from namelist parameters
312      !!----------------------------------------------------------------------
313      CHARACTER(len=6)        , INTENT(in   ) ::   cd_topbot       ! top/ bot indicator
314      INTEGER , DIMENSION(:,:), INTENT(in   ) ::   k_mk            ! 1st/last  wet level
315      REAL(wp)                , INTENT(  out) ::   pCdmin, pCdmax  ! min and max drag coef. [-]
316      REAL(wp)                , INTENT(  out) ::   pz0             ! roughness              [m]
317! tipaccs 2d top tidal velocity
318      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pke0            ! background KE          [m2/s2]
319! end tipaccs 2d top tidal velocity
320      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCd0            ! masked precomputed part of the non-linear drag coefficient
321      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCdU            ! minus linear drag*|U| at t-points  [m/s]
322      !!
323      CHARACTER(len=40) ::   cl_namdrg, cl_file, cl_varname, cl_namref, cl_namcfg  ! local names
324      INTEGER ::   ji, jj              ! dummy loop indexes
325      LOGICAL ::   ll_top, ll_bot      ! local logical
326      INTEGER ::   ios, inum, imk      ! local integers
327      REAL(wp)::   zmsk, zzz, zcd      ! local scalars
328      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk_boost   ! 2D workspace
329      !!
330      NAMELIST/namdrg_top/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost
331      NAMELIST/namdrg_bot/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost
332
333! tipaccs 2d top tidal velocity
334      CHARACTER(len=255) :: cn_dirttv
335      TYPE(FLD_N)        :: sn_ttv
336      REAL(wp)           :: rn_ttv_sf, rn_ttv_os
337      NAMELIST/namdrg_top_tipaccs/ ln_2d_ttv, rn_ttv_sf, rn_ttv_os, cn_dirttv, sn_ttv
338! end tipaccs 2d top tidal velocity
339      !!----------------------------------------------------------------------
340      !
341      !                          !==  set TOP / BOTTOM specificities  ==!
342      ll_top = .FALSE.
343      ll_bot = .FALSE.
344      !
345      SELECT CASE (cd_topbot)
346      CASE( 'TOP   ' )
347         ll_top = .TRUE.
348         cl_namdrg  = 'namdrg_top'
349         cl_namref  = 'namdrg_top in reference     namelist'
350         cl_namcfg  = 'namdrg_top in configuration namelist'
351         cl_file    = 'tfr_coef.nc'
352         cl_varname = 'tfr_coef'
353      CASE( 'BOTTOM' )
354         ll_bot = .TRUE.
355         cl_namdrg  = 'namdrg_bot'
356         cl_namref  = 'namdrg_bot  in reference     namelist'
357         cl_namcfg  = 'namdrg_bot  in configuration namelist'
358         cl_file    = 'bfr_coef.nc'
359         cl_varname = 'bfr_coef'
360      CASE DEFAULT
361         CALL ctl_stop( 'drg_init: bad value for cd_topbot ' )
362      END SELECT
363      !
364      !                          !==  read namlist  ==!
365      !
366      IF(ll_top)   READ  ( numnam_ref, namdrg_top, IOSTAT = ios, ERR = 901)
367      IF(ll_bot)   READ  ( numnam_ref, namdrg_bot, IOSTAT = ios, ERR = 901)
368901   IF( ios /= 0 )   CALL ctl_nam( ios , TRIM(cl_namref) )
369      IF(ll_top)   READ  ( numnam_cfg, namdrg_top, IOSTAT = ios, ERR = 902 )
370      IF(ll_bot)   READ  ( numnam_cfg, namdrg_bot, IOSTAT = ios, ERR = 902 )
371902   IF( ios >  0 )   CALL ctl_nam( ios , TRIM(cl_namcfg) )
372      IF(lwm .AND. ll_top)   WRITE ( numond, namdrg_top )
373      IF(lwm .AND. ll_bot)   WRITE ( numond, namdrg_bot )
374
375! tipaccs 2d top tidal velocity
376      IF (ll_top) THEN
377         READ (numnam_ref, namdrg_top_tipaccs, IOSTAT = ios, ERR = 905)
378905      IF( ios /= 0 )   CALL ctl_nam( ios , TRIM(cl_namref)//'_tipaccs' )
379         READ (numnam_cfg, namdrg_top_tipaccs, IOSTAT = ios, ERR = 906)
380906      IF( ios /= 0 )   CALL ctl_nam( ios , TRIM(cl_namcfg)//'_tipaccs' )
381      ENDIF
382! end tipaccs 2d top tidal velocity
383
384      !
385      IF(lwp) THEN
386         WRITE(numout,*)
387         WRITE(numout,*) '   Namelist ',TRIM(cl_namdrg),' : set ',TRIM(cd_topbot),' friction parameters'
388         WRITE(numout,*) '      drag coefficient                        rn_Cd0   = ', rn_Cd0
389         WRITE(numout,*) '      characteristic velocity (linear case)   rn_Uc0   = ', rn_Uc0, ' m/s'
390         WRITE(numout,*) '      non-linear drag maximum                 rn_Cdmax = ', rn_Cdmax
391
392! tipaccs 2d top tidal velocity
393         IF (ll_top) WRITE(numout,*) '      use 2d top tidal velocity               ln_2d_ttv= ', ln_2d_ttv
394         IF (ln_2d_ttv .AND. ll_top) THEN
395            WRITE(numout,*) '      2d top tital velocity read from file    sn_ttv   = ',TRIM(sn_ttv%clname)
396            WRITE(numout,*) '      scale factor applied to ttv             rn_ttv_sf= ', rn_ttv_sf
397            WRITE(numout,*) '      offset       applied to ttv             rn_ttv_os= ', rn_ttv_os
398         ELSE
399            WRITE(numout,*) '      background kinetic energy  (n-l case)   rn_ke0   = ', rn_ke0
400         ENDIF
401! end tipaccs 2d top tidal velocity
402         WRITE(numout,*) '      bottom roughness           (n-l case)   rn_z0    = ', rn_z0
403         WRITE(numout,*) '      set a regional boost of Cd0             ln_boost = ', ln_boost
404         WRITE(numout,*) '         associated boost factor              rn_boost = ', rn_boost
405      ENDIF
406      !
407      !                          !==  return some namelist parametres  ==!   (used in non_lin and loglayer cases)
408      pCdmin = rn_Cd0
409      pCdmax = rn_Cdmax
410      pz0    = rn_z0
411      pke0   = rn_ke0
412      !
413! tipaccs 2d top tidal velocity
414      IF (ln_2d_ttv) THEN
415                  IF(lwp) WRITE(numout,*)
416         IF(lwp) WRITE(numout,*) '   ==>>>   use a 2d top tidal velocity read in ',TRIM(sn_ttv%clname), ' file'
417         ! 2d top tidal velocity
418         CALL iom_open ( TRIM(sn_ttv%clname), inum )
419         CALL iom_get  ( inum, jpdom_global, TRIM(sn_ttv%clvar), pke0, 1 )
420         CALL iom_close( inum)
421         !
422         ! Eq 9c Jourdain et al. (2019)
423         ! input file is a velocity, NEMO need a square velocity
424         pke0 = pke0 * rn_ttv_sf + rn_ttv_os
425         pke0 = pke0 * pke0
426         !
427      ENDIF
428! end tipaccs 2d top tidal velocity
429      !                          !==  mask * boost factor  ==!
430      !
431      IF( ln_boost ) THEN           !* regional boost:   boost factor = 1 + regional boost
432         IF(lwp) WRITE(numout,*)
433         IF(lwp) WRITE(numout,*) '   ==>>>   use a regional boost read in ', TRIM(cl_file), ' file'
434         IF(lwp) WRITE(numout,*) '           using enhancement factor of ', rn_boost
435         ! cl_varname is a coefficient in [0,1] giving where to apply the regional boost
436         CALL iom_open ( TRIM(cl_file), inum )
437         CALL iom_get  ( inum, jpdom_global, TRIM(cl_varname), zmsk_boost, 1 )
438         CALL iom_close( inum)
439         zmsk_boost(:,:) = 1._wp + rn_boost * zmsk_boost(:,:)
440         !
441      ELSE                          !* no boost:   boost factor = 1
442         zmsk_boost(:,:) = 1._wp
443      ENDIF
444      !                             !* mask outside ocean cavities area (top) or land area (bot)
445      IF(ll_top)   zmsk_boost(:,:) = zmsk_boost(:,:) * ssmask(:,:) * (1. - tmask(:,:,1) )  ! none zero in ocean cavities only
446      IF(ll_bot)   zmsk_boost(:,:) = zmsk_boost(:,:) * ssmask(:,:)                         ! x seafloor mask
447      !
448      l_log_not_linssh = .FALSE.    ! default definition
449      !
450      SELECT CASE( ndrg )
451      !
452      CASE( np_OFF  )            !==  No top/bottom friction  ==!   (pCdU = 0)
453         IF(lwp) WRITE(numout,*)
454         IF(lwp) WRITE(numout,*) '   ==>>>   ',TRIM(cd_topbot),' free-slip, friction set to zero'
455         !
456         l_zdfdrg = .FALSE.         ! no time variation of the drag: set it one for all
457         !
458         pCdU(:,:) = 0._wp
459         pCd0(:,:) = 0._wp
460         !
461      CASE( np_lin )             !==  linear friction  ==!   (pCdU = Cd0 * Uc0)
462         IF(lwp) WRITE(numout,*)
463         IF(lwp) WRITE(numout,*) '   ==>>>   linear ',TRIM(cd_topbot),' friction (constant coef = Cd0*Uc0 = ', rn_Cd0*rn_Uc0, ')'
464         !
465         l_zdfdrg = .FALSE.         ! no time variation of the Cd*|U| : set it one for all
466         !                     
467         pCd0(:,:) = rn_Cd0 * zmsk_boost(:,:)  !* constant in time drag coefficient (= mask (and boost) Cd0)
468         pCdU(:,:) = - pCd0(:,:) * rn_Uc0      !  using a constant velocity
469         !
470      CASE( np_non_lin )         !== non-linear friction  ==!   (pCd0 = Cd0 )
471         IF(lwp) WRITE(numout,*)
472         IF(lwp) WRITE(numout,*) '   ==>>>   quadratic ',TRIM(cd_topbot),' friction (propotional to module of the velocity)'
473         IF(lwp) WRITE(numout,*) '   with    a drag coefficient Cd0 = ', rn_Cd0, ', and'
474! tipaccs 2d top tidal velocity
475         IF (ln_2d_ttv) THEN
476            IF(lwp) WRITE(numout,*) '           a 2d top tidal velocity from',TRIM(sn_ttv%clname)
477         ELSE
478            IF(lwp) WRITE(numout,*) '           a background velocity module of (rn_ke0)^1/2 = ', SQRT(rn_ke0), 'm/s)'
479         END IF
480! end tipaccs 2d top tidal velocity
481         !
482         l_zdfdrg = .TRUE.          !* Cd*|U| updated at each time-step (it depends on ocean velocity)
483         !
484         pCd0(:,:) = rn_Cd0 * zmsk_boost(:,:)  !* constant in time proportionality coefficient (= mask (and boost) Cd0)
485         pCdU(:,:) = 0._wp                     
486         !
487      CASE( np_loglayer )       !== logarithmic layer formulation of friction  ==!   (CdU = (vkarman log(z/z0))^2 |U| )
488         IF(lwp) WRITE(numout,*)
489         IF(lwp) WRITE(numout,*) '   ==>>>   quadratic ',TRIM(cd_topbot),' drag (propotional to module of the velocity)'
490         IF(lwp) WRITE(numout,*) '   with   a logarithmic Cd0 formulation Cd0 = ( vkarman log(z/z0) )^2 ,'
491! tipaccs 2d top tidal velocity
492         IF (ln_2d_ttv) THEN
493            IF(lwp) WRITE(numout,*) '           a 2d top tidal velocity from',TRIM(sn_ttv%clname)
494         ELSE
495            IF(lwp) WRITE(numout,*) '           a background velocity module of (rn_ke0)^1/2 = ', SQRT(rn_ke0), 'm/s)'
496         END IF
497! tipaccs 2d top tidal velocity
498         IF(lwp) WRITE(numout,*) '          a logarithmic formulation: a roughness of ', pz0, ' meters,   and '
499         IF(lwp) WRITE(numout,*) '          a proportionality factor bounded by min/max values of ', pCdmin, pCdmax
500         !
501         l_zdfdrg = .TRUE.          !* Cd*|U| updated at each time-step (it depends on ocean velocity)
502         !
503         IF( ln_linssh ) THEN       !* pCd0 = (v log(z/z0))^2   as velocity points have a fixed z position
504            IF(lwp) WRITE(numout,*)
505            IF(lwp) WRITE(numout,*) '   N.B.   linear free surface case, Cd0 computed one for all'
506            !
507            l_log_not_linssh = .FALSE.    !- don't update Cd at each time step
508            !
509            DO_2D( 1, 1, 1, 1 )              ! pCd0 = mask (and boosted) logarithmic drag coef.
510               zzz =  0.5_wp * e3t_0(ji,jj,k_mk(ji,jj))
511               zcd = (  vkarmn / LOG( zzz / rn_z0 )  )**2
512               pCd0(ji,jj) = zmsk_boost(ji,jj) * MIN(  MAX( rn_Cd0 , zcd ) , rn_Cdmax  )  ! rn_Cd0 < Cd0 < rn_Cdmax
513            END_2D
514         ELSE                       !* Cd updated at each time-step ==> pCd0 = mask * boost
515            IF(lwp) WRITE(numout,*)
516            IF(lwp) WRITE(numout,*) '   N.B.   non-linear free surface case, Cd0 updated at each time-step '
517            !
518            l_log_not_linssh = .TRUE.     ! compute the drag coef. at each time-step
519            !
520            pCd0(:,:) = zmsk_boost(:,:)
521         ENDIF
522         pCdU(:,:) = 0._wp          ! initialisation to zero (will be updated at each time step)
523         !
524      CASE DEFAULT
525         CALL ctl_stop( 'drg_init: bad flag value for ndrg ' )
526      END SELECT
527      !
528   END SUBROUTINE drg_init
529
530   !!======================================================================
531END MODULE zdfdrg
Note: See TracBrowser for help on using the repository browser.