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, 16 months 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.