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

source: NEMO/trunk/src/OCE/ZDF/zdfdrg.F90 @ 12377

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

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 24.4 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 "do_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, Kmm, 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      INTEGER                 , INTENT(in   ) ::   Kmm      ! ocean time level index
102      !                       !               !!         !==  top or bottom variables  ==!
103      INTEGER , DIMENSION(:,:), INTENT(in   ) ::   k_mk     ! wet level (1st or last)
104      REAL(wp)                , INTENT(in   ) ::   pCdmin   ! min drag value
105      REAL(wp)                , INTENT(in   ) ::   pCdmax   ! max drag value
106      REAL(wp)                , INTENT(in   ) ::   pz0      ! roughness
107      REAL(wp)                , INTENT(in   ) ::   pke0     ! background tidal KE
108      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pCd0     ! masked precomputed part of Cd0
109      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCdU     ! = - Cd*|U|   (t-points) [m/s]
110      !!
111      INTEGER ::   ji, jj   ! dummy loop indices
112      INTEGER ::   imk      ! local integers
113      REAL(wp)::   zzz, zut, zvt, zcd   ! local scalars
114      !!----------------------------------------------------------------------
115      !
116      IF( l_log_not_linssh ) THEN     !==  "log layer"  ==!   compute Cd and -Cd*|U|
117         DO_2D_00_00
118            imk = k_mk(ji,jj)          ! ocean bottom level at t-points
119            zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point
120            zvt = vv(ji,jj,imk,Kmm) + vv(ji,jj-1,imk,Kmm)
121            zzz = 0.5_wp * e3t(ji,jj,imk,Kmm)           ! 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_2D
128      ELSE                                            !==  standard Cd  ==!
129         DO_2D_00_00
130            imk = k_mk(ji,jj)    ! ocean bottom level at t-points
131            zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point
132            zvt = vv(ji,jj,imk,Kmm) + vv(ji,jj-1,imk,Kmm)
133            !                                                           ! here pCd0 = mask*boost * drag
134            pCdU(ji,jj) = - pCd0(ji,jj) * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  )
135         END_2D
136      ENDIF
137      !
138      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pCdU, clinfo1=' Cd*U ')
139      !
140   END SUBROUTINE zdf_drg
141
142
143   SUBROUTINE zdf_drg_exp( kt, Kmm, pub, pvb, pua, pva )
144      !!----------------------------------------------------------------------
145      !!                  ***  ROUTINE zdf_drg_exp  ***
146      !!
147      !! ** Purpose :   compute and add the explicit top and bottom frictions.
148      !!
149      !! ** Method  :   in explicit case,
150      !!
151      !!              NB: in implicit case the calculation is performed in dynzdf.F90
152      !!
153      !! ** Action  :   (pua,pva)   momentum trend increased by top & bottom friction trend
154      !!---------------------------------------------------------------------
155      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
156      INTEGER                         , INTENT(in   ) ::   Kmm        ! time level indices
157      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pub, pvb   ! the two components of the before velocity
158      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! the two components of the velocity tendency
159      !!
160      INTEGER  ::   ji, jj       ! dummy loop indexes
161      INTEGER  ::   ikbu, ikbv   ! local integers
162      REAL(wp) ::   zm1_2dt      ! local scalar
163      REAL(wp) ::   zCdu, zCdv   !   -      -
164      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv
165      !!---------------------------------------------------------------------
166      !
167!!gm bug : time step is only rdt (not 2 rdt if euler start !)
168      zm1_2dt = - 1._wp / ( 2._wp * rdt )
169
170      IF( l_trddyn ) THEN      ! trends: store the input trends
171         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
172         ztrdu(:,:,:) = pua(:,:,:)
173         ztrdv(:,:,:) = pva(:,:,:)
174      ENDIF
175
176      DO_2D_00_00
177         ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
178         ikbv = mbkv(ji,jj)
179         !
180         ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
181         zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u(ji,jj,ikbu,Kmm)
182         zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v(ji,jj,ikbv,Kmm)
183         !
184         pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
185         pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv)
186      END_2D
187      !
188      IF( ln_isfcav ) THEN        ! ocean cavities
189         DO_2D_00_00
190            ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
191            ikbv = mikv(ji,jj)
192            !
193            ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
194            zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u(ji,jj,ikbu,Kmm)    ! NB: Cdtop masked
195            zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v(ji,jj,ikbv,Kmm)
196            !
197            pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu)
198            pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv)
199         END_2D
200      ENDIF
201      !
202      IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics
203         ztrdu(:,:,:) = pua(:,:,:) - ztrdu(:,:,:)
204         ztrdv(:,:,:) = pva(:,:,:) - ztrdv(:,:,:)
205         CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt, Kmm )
206         DEALLOCATE( ztrdu, ztrdv )
207      ENDIF
208      !                                          ! print mean trends (used for debugging)
209      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pua, clinfo1=' bfr  - Ua: ', mask1=umask,               &
210         &                                  tab3d_2=pva, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
211      !
212   END SUBROUTINE zdf_drg_exp
213
214
215   SUBROUTINE zdf_drg_init
216      !!----------------------------------------------------------------------
217      !!                  ***  ROUTINE zdf_brg_init  ***
218      !!
219      !! ** Purpose :   Initialization of the bottom friction
220      !!
221      !! ** Method  :   Read the namdrg namelist and check their consistency
222      !!                called at the first timestep (nit000)
223      !!----------------------------------------------------------------------
224      INTEGER   ::   ji, jj      ! dummy loop indexes
225      INTEGER   ::   ios, ioptio   ! local integers
226      !!
227      NAMELIST/namdrg/ ln_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp
228      !!----------------------------------------------------------------------
229      !
230      !                     !==  drag nature  ==!
231      !
232      READ  ( numnam_ref, namdrg, IOSTAT = ios, ERR = 901)
233901   IF( ios /= 0 )   CALL ctl_nam( ios , 'namdrg in reference namelist' )
234      READ  ( numnam_cfg, namdrg, IOSTAT = ios, ERR = 902 )
235902   IF( ios >  0 )   CALL ctl_nam( ios , 'namdrg in configuration namelist' )
236      IF(lwm) WRITE ( numond, namdrg )
237      !
238      IF(lwp) THEN
239         WRITE(numout,*)
240         WRITE(numout,*) 'zdf_drg_init : top and/or bottom drag setting'
241         WRITE(numout,*) '~~~~~~~~~~~~'
242         WRITE(numout,*) '   Namelist namdrg : top/bottom friction choices'
243         WRITE(numout,*) '      free-slip       : Cd = 0                  ln_OFF      = ', ln_OFF 
244         WRITE(numout,*) '      linear  drag    : Cd = Cd0                ln_lin      = ', ln_lin
245         WRITE(numout,*) '      non-linear  drag: Cd = Cd0_nl |U|         ln_non_lin  = ', ln_non_lin
246         WRITE(numout,*) '      logarithmic drag: Cd = vkarmn/log(z/z0)   ln_loglayer = ', ln_loglayer
247         WRITE(numout,*) '      implicit friction                         ln_drgimp   = ', ln_drgimp
248      ENDIF
249      !
250      ioptio = 0                       ! set ndrg and control check
251      IF( ln_OFF      ) THEN   ;   ndrg = np_OFF        ;   ioptio = ioptio + 1   ;   ENDIF
252      IF( ln_lin      ) THEN   ;   ndrg = np_lin        ;   ioptio = ioptio + 1   ;   ENDIF
253      IF( ln_non_lin  ) THEN   ;   ndrg = np_non_lin    ;   ioptio = ioptio + 1   ;   ENDIF
254      IF( ln_loglayer ) THEN   ;   ndrg = np_loglayer   ;   ioptio = ioptio + 1   ;   ENDIF
255      !
256      IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' )
257      !
258      !
259      !                     !==  BOTTOM drag setting  ==!   (applied at seafloor)
260      !
261      ALLOCATE( rCd0_bot(jpi,jpj), rCdU_bot(jpi,jpj) )
262      CALL drg_init( 'BOTTOM'   , mbkt       ,                                         &   ! <== in
263         &           r_Cdmin_bot, r_Cdmax_bot, r_z0_bot, r_ke0_bot, rCd0_bot, rCdU_bot )   ! ==> out
264
265      !
266      !                     !==  TOP drag setting  ==!   (applied at the top of ocean cavities)
267      !
268      IF( ln_isfcav ) THEN              ! Ocean cavities: top friction setting
269         ALLOCATE( rCd0_top(jpi,jpj), rCdU_top(jpi,jpj) )
270         CALL drg_init( 'TOP   '   , mikt       ,                                         &   ! <== in
271            &           r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top, rCd0_top, rCdU_top )   ! ==> out
272      ENDIF
273      !
274   END SUBROUTINE zdf_drg_init
275
276
277   SUBROUTINE drg_init( cd_topbot, k_mk,  &
278      &                 pCdmin, pCdmax, pz0, pke0, pCd0, pCdU ) 
279      !!----------------------------------------------------------------------
280      !!                  ***  ROUTINE drg_init  ***
281      !!
282      !! ** Purpose :   Initialization of the top/bottom friction CdO and Cd
283      !!              from namelist parameters
284      !!----------------------------------------------------------------------
285      CHARACTER(len=6)        , INTENT(in   ) ::   cd_topbot       ! top/ bot indicator
286      INTEGER , DIMENSION(:,:), INTENT(in   ) ::   k_mk            ! 1st/last  wet level
287      REAL(wp)                , INTENT(  out) ::   pCdmin, pCdmax  ! min and max drag coef. [-]
288      REAL(wp)                , INTENT(  out) ::   pz0             ! roughness              [m]
289      REAL(wp)                , INTENT(  out) ::   pke0            ! background KE          [m2/s2]
290      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCd0            ! masked precomputed part of the non-linear drag coefficient
291      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCdU            ! minus linear drag*|U| at t-points  [m/s]
292      !!
293      CHARACTER(len=40) ::   cl_namdrg, cl_file, cl_varname, cl_namref, cl_namcfg  ! local names
294      INTEGER ::   ji, jj              ! dummy loop indexes
295      LOGICAL ::   ll_top, ll_bot      ! local logical
296      INTEGER ::   ios, inum, imk      ! local integers
297      REAL(wp)::   zmsk, zzz, zcd      ! local scalars
298      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk_boost   ! 2D workspace
299      !!
300      NAMELIST/namdrg_top/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost
301      NAMELIST/namdrg_bot/ rn_Cd0, rn_Uc0, rn_Cdmax, rn_ke0, rn_z0, ln_boost, rn_boost
302      !!----------------------------------------------------------------------
303      !
304      !                          !==  set TOP / BOTTOM specificities  ==!
305      ll_top = .FALSE.
306      ll_bot = .FALSE.
307      !
308      SELECT CASE (cd_topbot)
309      CASE( 'TOP   ' )
310         ll_top = .TRUE.
311         cl_namdrg  = 'namdrg_top'
312         cl_namref  = 'namdrg_top in reference     namelist'
313         cl_namcfg  = 'namdrg_top in configuration namelist'
314         cl_file    = 'tfr_coef.nc'
315         cl_varname = 'tfr_coef'
316      CASE( 'BOTTOM' )
317         ll_bot = .TRUE.
318         cl_namdrg  = 'namdrg_bot'
319         cl_namref  = 'namdrg_bot  in reference     namelist'
320         cl_namcfg  = 'namdrg_bot  in configuration namelist'
321         cl_file    = 'bfr_coef.nc'
322         cl_varname = 'bfr_coef'
323      CASE DEFAULT
324         CALL ctl_stop( 'drg_init: bad value for cd_topbot ' )
325      END SELECT
326      !
327      !                          !==  read namlist  ==!
328      !
329      IF(ll_top)   READ  ( numnam_ref, namdrg_top, IOSTAT = ios, ERR = 901)
330      IF(ll_bot)   READ  ( numnam_ref, namdrg_bot, IOSTAT = ios, ERR = 901)
331901   IF( ios /= 0 )   CALL ctl_nam( ios , TRIM(cl_namref) )
332      IF(ll_top)   READ  ( numnam_cfg, namdrg_top, IOSTAT = ios, ERR = 902 )
333      IF(ll_bot)   READ  ( numnam_cfg, namdrg_bot, IOSTAT = ios, ERR = 902 )
334902   IF( ios >  0 )   CALL ctl_nam( ios , TRIM(cl_namcfg) )
335      IF(lwm .AND. ll_top)   WRITE ( numond, namdrg_top )
336      IF(lwm .AND. ll_bot)   WRITE ( numond, namdrg_bot )
337      !
338      IF(lwp) THEN
339         WRITE(numout,*)
340         WRITE(numout,*) '   Namelist ',TRIM(cl_namdrg),' : set ',TRIM(cd_topbot),' friction parameters'
341         WRITE(numout,*) '      drag coefficient                        rn_Cd0   = ', rn_Cd0
342         WRITE(numout,*) '      characteristic velocity (linear case)   rn_Uc0   = ', rn_Uc0, ' m/s'
343         WRITE(numout,*) '      non-linear drag maximum                 rn_Cdmax = ', rn_Cdmax
344         WRITE(numout,*) '      background kinetic energy  (n-l case)   rn_ke0   = ', rn_ke0
345         WRITE(numout,*) '      bottom roughness           (n-l case)   rn_z0    = ', rn_z0
346         WRITE(numout,*) '      set a regional boost of Cd0             ln_boost = ', ln_boost
347         WRITE(numout,*) '         associated boost factor              rn_boost = ', rn_boost
348      ENDIF
349      !
350      !                          !==  return some namelist parametres  ==!   (used in non_lin and loglayer cases)
351      pCdmin = rn_Cd0
352      pCdmax = rn_Cdmax
353      pz0    = rn_z0
354      pke0   = rn_ke0
355      !
356      !                          !==  mask * boost factor  ==!
357      !
358      IF( ln_boost ) THEN           !* regional boost:   boost factor = 1 + regional boost
359         IF(lwp) WRITE(numout,*)
360         IF(lwp) WRITE(numout,*) '   ==>>>   use a regional boost read in ', TRIM(cl_file), ' file'
361         IF(lwp) WRITE(numout,*) '           using enhancement factor of ', rn_boost
362         ! cl_varname is a coefficient in [0,1] giving where to apply the regional boost
363         CALL iom_open ( TRIM(cl_file), inum )
364         CALL iom_get  ( inum, jpdom_data, TRIM(cl_varname), zmsk_boost, 1 )
365         CALL iom_close( inum)
366         zmsk_boost(:,:) = 1._wp + rn_boost * zmsk_boost(:,:)
367         !
368      ELSE                          !* no boost:   boost factor = 1
369         zmsk_boost(:,:) = 1._wp
370      ENDIF
371      !                             !* mask outside ocean cavities area (top) or land area (bot)
372      IF(ll_top)   zmsk_boost(:,:) = zmsk_boost(:,:) * ssmask(:,:) * (1. - tmask(:,:,1) )  ! none zero in ocean cavities only
373      IF(ll_bot)   zmsk_boost(:,:) = zmsk_boost(:,:) * ssmask(:,:)                         ! x seafloor mask
374      !
375      !
376      SELECT CASE( ndrg )
377      !
378      CASE( np_OFF  )            !==  No top/bottom friction  ==!   (pCdU = 0)
379         IF(lwp) WRITE(numout,*)
380         IF(lwp) WRITE(numout,*) '   ==>>>   ',TRIM(cd_topbot),' free-slip, friction set to zero'
381         !
382         l_zdfdrg = .FALSE.         ! no time variation of the drag: set it one for all
383         !
384         pCdU(:,:) = 0._wp
385         pCd0(:,:) = 0._wp
386         !
387      CASE( np_lin )             !==  linear friction  ==!   (pCdU = Cd0 * Uc0)
388         IF(lwp) WRITE(numout,*)
389         IF(lwp) WRITE(numout,*) '   ==>>>   linear ',TRIM(cd_topbot),' friction (constant coef = Cd0*Uc0 = ', rn_Cd0*rn_Uc0, ')'
390         !
391         l_zdfdrg = .FALSE.         ! no time variation of the Cd*|U| : set it one for all
392         !                     
393         pCd0(:,:) = rn_Cd0 * zmsk_boost(:,:)  !* constant in time drag coefficient (= mask (and boost) Cd0)
394         pCdU(:,:) = - pCd0(:,:) * rn_Uc0      !  using a constant velocity
395         !
396      CASE( np_non_lin )         !== non-linear friction  ==!   (pCd0 = Cd0 )
397         IF(lwp) WRITE(numout,*)
398         IF(lwp) WRITE(numout,*) '   ==>>>   quadratic ',TRIM(cd_topbot),' friction (propotional to module of the velocity)'
399         IF(lwp) WRITE(numout,*) '   with    a drag coefficient Cd0 = ', rn_Cd0, ', and'
400         IF(lwp) WRITE(numout,*) '           a background velocity module of (rn_ke0)^1/2 = ', SQRT(rn_ke0), 'm/s)'
401         !
402         l_zdfdrg = .TRUE.          !* Cd*|U| updated at each time-step (it depends on ocean velocity)
403         !
404         pCd0(:,:) = rn_Cd0 * zmsk_boost(:,:)  !* constant in time proportionality coefficient (= mask (and boost) Cd0)
405         pCdU(:,:) = 0._wp                     ! 
406         !
407      CASE( np_loglayer )       !== logarithmic layer formulation of friction  ==!   (CdU = (vkarman log(z/z0))^2 |U| )
408         IF(lwp) WRITE(numout,*)
409         IF(lwp) WRITE(numout,*) '   ==>>>   quadratic ',TRIM(cd_topbot),' drag (propotional to module of the velocity)'
410         IF(lwp) WRITE(numout,*) '   with   a logarithmic Cd0 formulation Cd0 = ( vkarman log(z/z0) )^2 ,'
411         IF(lwp) WRITE(numout,*) '          a background velocity module of (rn_ke0)^1/2 = ', SQRT(pke0), 'm/s), '
412         IF(lwp) WRITE(numout,*) '          a logarithmic formulation: a roughness of ', pz0, ' meters,   and '
413         IF(lwp) WRITE(numout,*) '          a proportionality factor bounded by min/max values of ', pCdmin, pCdmax
414         !
415         l_zdfdrg = .TRUE.          !* Cd*|U| updated at each time-step (it depends on ocean velocity)
416         !
417         IF( ln_linssh ) THEN       !* pCd0 = (v log(z/z0))^2   as velocity points have a fixed z position
418            IF(lwp) WRITE(numout,*)
419            IF(lwp) WRITE(numout,*) '   N.B.   linear free surface case, Cd0 computed one for all'
420            !
421            l_log_not_linssh = .FALSE.    !- don't update Cd at each time step
422            !
423            DO_2D_11_11
424               zzz =  0.5_wp * e3t_0(ji,jj,k_mk(ji,jj))
425               zcd = (  vkarmn / LOG( zzz / rn_z0 )  )**2
426               pCd0(ji,jj) = zmsk_boost(ji,jj) * MIN(  MAX( rn_Cd0 , zcd ) , rn_Cdmax  )  ! rn_Cd0 < Cd0 < rn_Cdmax
427            END_2D
428         ELSE                       !* Cd updated at each time-step ==> pCd0 = mask * boost
429            IF(lwp) WRITE(numout,*)
430            IF(lwp) WRITE(numout,*) '   N.B.   non-linear free surface case, Cd0 updated at each time-step '
431            !
432            l_log_not_linssh = .TRUE.     ! compute the drag coef. at each time-step
433            !
434            pCd0(:,:) = zmsk_boost(:,:)
435         ENDIF
436         pCdU(:,:) = 0._wp          ! initialisation to zero (will be updated at each time step)
437         !
438      CASE DEFAULT
439         CALL ctl_stop( 'drg_init: bad flag value for ndrg ' )
440      END SELECT
441      !
442   END SUBROUTINE drg_init
443
444   !!======================================================================
445END MODULE zdfdrg
Note: See TracBrowser for help on using the repository browser.