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.
zdfbfr.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90 @ 4821

Last change on this file since 4821 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
  • Property svn:keywords set to Id
File size: 16.6 KB
Line 
1MODULE zdfbfr
2   !!======================================================================
3   !!                       ***  MODULE  zdfbfr  ***
4   !! Ocean physics: 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   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   zdf_bfr      : update bottom friction coefficient (non-linear bottom friction only)
16   !!   zdf_bfr_init : read in namelist and control the bottom friction parameters.
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE zdf_oce         ! ocean vertical physics variables
21   USE in_out_manager  ! I/O manager
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE lib_mpp         ! distributed memory computing
24   USE prtctl          ! Print control
25   USE timing          ! Timing
26   USE wrk_nemo        ! Memory Allocation
27   USE phycst, ONLY: vkarmn
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   zdf_bfr         ! called by step.F90
33   PUBLIC   zdf_bfr_init    ! called by nemogcm.F90
34
35   !                                 !!* Namelist nambfr: bottom friction namelist *
36   INTEGER , PUBLIC ::   nn_bfr       ! = 0/1/2/3 type of bottom friction  (PUBLIC for TAM)
37   REAL(wp), PUBLIC ::   rn_bfri1     ! bottom drag coefficient (linear case)  (PUBLIC for TAM)
38   REAL(wp), PUBLIC ::   rn_bfri2     ! bottom drag coefficient (non linear case) (PUBLIC for TAM)
39   REAL(wp), PUBLIC ::   rn_bfri2_max ! Maximum bottom drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM)
40   REAL(wp), PUBLIC ::   rn_bfeb2     ! background bottom turbulent kinetic energy  [m2/s2] (PUBLIC for TAM)
41   REAL(wp), PUBLIC ::   rn_bfrien    ! local factor to enhance coefficient bfri (PUBLIC for TAM)
42   LOGICAL , PUBLIC ::   ln_bfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM)
43   LOGICAL , PUBLIC ::   ln_loglayer  ! switch for log layer bfr coeff. (PUBLIC for TAM)
44   REAL(wp), PUBLIC ::   rn_bfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
45   LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d            ! 2D bottom drag coefficient (PUBLIC for TAM)
47
48   !! * Substitutions
49#  include "vectopt_loop_substitute.h90"
50#  include "domzgr_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
53   !! $Id$
54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   INTEGER FUNCTION zdf_bfr_alloc()
59      !!----------------------------------------------------------------------
60      !!                ***  FUNCTION zdf_bfr_alloc  ***
61      !!----------------------------------------------------------------------
62      ALLOCATE( bfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
63      !
64      IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc )
65      IF( zdf_bfr_alloc /= 0 )   CALL ctl_warn('zdf_bfr_alloc: failed to allocate arrays.')
66   END FUNCTION zdf_bfr_alloc
67
68
69   SUBROUTINE zdf_bfr( kt )
70      !!----------------------------------------------------------------------
71      !!                   ***  ROUTINE zdf_bfr  ***
72      !!
73      !! ** Purpose :   compute the bottom friction coefficient.
74      !!
75      !! ** Method  :   Calculate and store part of the momentum trend due
76      !!              to bottom friction following the chosen friction type
77      !!              (free-slip, linear, or quadratic). The component
78      !!              calculated here is multiplied by the bottom velocity in
79      !!              dyn_bfr to provide the trend term.
80      !!                The coefficients are updated at each time step only
81      !!              in the quadratic case.
82      !!
83      !! ** Action  :   bfrua , bfrva   bottom friction coefficients
84      !!----------------------------------------------------------------------
85      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
86      !!
87      INTEGER  ::   ji, jj                       ! dummy loop indices
88      INTEGER  ::   ikbt, ikbu, ikbv             ! local integers
89      REAL(wp) ::   zvu, zuv, zecu, zecv, ztmp   ! temporary scalars
90      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt
91      !!----------------------------------------------------------------------
92      !
93      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr')
94      !
95      IF( kt == nit000 .AND. lwp ) THEN
96         WRITE(numout,*)
97         WRITE(numout,*) 'zdf_bfr : Set bottom friction coefficient (non-linear case)'
98         WRITE(numout,*) '~~~~~~~~'
99      ENDIF
100      !
101      IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only
102         !
103         CALL wrk_alloc( jpi, jpj, zbfrt )
104
105         IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient
106
107#  if defined key_vectopt_loop
108            DO jj = 1, 1
109!CDIR NOVERRCHK
110               DO ji = 1, jpij   ! vector opt. (forced unrolling)
111#  else
112!CDIR NOVERRCHK
113            DO jj = 1, jpj
114!CDIR NOVERRCHK
115               DO ji = 1, jpi
116#  endif
117                  ikbt = mbkt(ji,jj)
118! JC: possible WAD implementation should modify line below if layers vanish
119                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
120                  zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
121                  zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max)
122               END DO
123            END DO
124         !   
125         ELSE
126            zbfrt(:,:) = bfrcoef2d(:,:)
127         ENDIF
128
129# if defined key_vectopt_loop
130         DO jj = 1, 1
131!CDIR NOVERRCHK
132            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
133# else
134!CDIR NOVERRCHK
135         DO jj = 2, jpjm1
136!CDIR NOVERRCHK
137            DO ji = 2, jpim1
138# endif
139               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
140               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
141               !
142               zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
143                  &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
144               zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
145                  &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
146               !
147               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2  )
148               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  )
149               !
150               bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) * zecu
151               bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv
152            END DO
153         END DO
154
155         !
156         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
157         !
158         IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
159            &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 )
160         CALL wrk_dealloc( jpi,jpj,zbfrt )
161      ENDIF
162      !
163      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr')
164      !
165   END SUBROUTINE zdf_bfr
166
167
168   SUBROUTINE zdf_bfr_init
169      !!----------------------------------------------------------------------
170      !!                  ***  ROUTINE zdf_bfr_init  ***
171      !!
172      !! ** Purpose :   Initialization of the bottom friction
173      !!
174      !! ** Method  :   Read the nambfr namelist and check their consistency
175      !!                called at the first timestep (nit000)
176      !!----------------------------------------------------------------------
177      USE iom   ! I/O module for ehanced bottom friction file
178      !!
179      INTEGER   ::   inum         ! logical unit for enhanced bottom friction file
180      INTEGER   ::   ji, jj       ! dummy loop indexes
181      INTEGER   ::   ikbt, ikbu, ikbv   ! temporary integers
182      INTEGER   ::   ictu, ictv         !    -          -
183      INTEGER   ::   ios
184      REAL(wp)  ::   zminbfr, zmaxbfr   ! temporary scalars
185      REAL(wp)  ::   ztmp, zfru, zfrv   !    -         -
186      !!
187      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, &
188                    &  rn_bfrien, ln_bfrimp, ln_loglayer
189      !!----------------------------------------------------------------------
190      !
191      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr_init')
192      !
193      !                              !* Allocate zdfbfr arrays
194      IF( zdf_bfr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' )
195      !
196      !                              !* Parameter control and print
197      !
198      REWIND( numnam_ref )              ! Namelist nambfr in reference namelist : Bottom momentum boundary condition
199      READ  ( numnam_ref, nambfr, IOSTAT = ios, ERR = 901)
200901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambfr in reference namelist', lwp )
201
202      REWIND( numnam_cfg )              ! Namelist nambfr in configuration namelist : Bottom momentum boundary condition
203      READ  ( numnam_cfg, nambfr, IOSTAT = ios, ERR = 902 )
204902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambfr in configuration namelist', lwp )
205      IF(lwm) WRITE ( numond, nambfr )
206      IF(lwp) WRITE(numout,*)
207      IF(lwp) WRITE(numout,*) 'zdf_bfr_init : momentum bottom friction'
208      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
209      IF(lwp) WRITE(numout,*) '   Namelist nam_bfr : set bottom friction parameters'
210      !
211      SELECT CASE (nn_bfr)
212      !
213      CASE( 0 )
214         IF(lwp) WRITE(numout,*) '      free-slip '
215         bfrua(:,:) = 0.e0
216         bfrva(:,:) = 0.e0
217         !
218      CASE( 1 )
219         IF(lwp) WRITE(numout,*) '      linear botton friction'
220         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri1  = ', rn_bfri1
221         IF( ln_bfr2d ) THEN
222            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d
223            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
224         ENDIF
225         !
226         IF(ln_bfr2d) THEN
227            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
228            CALL iom_open('bfr_coef.nc',inum)
229            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
230            CALL iom_close(inum)
231            bfrcoef2d(:,:) = rn_bfri1 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
232         ELSE
233            bfrcoef2d(:,:) = rn_bfri1  ! initialize bfrcoef2d to the namelist variable
234         ENDIF
235         !
236         bfrua(:,:) = - bfrcoef2d(:,:)
237         bfrva(:,:) = - bfrcoef2d(:,:)
238         !
239      CASE( 2 )
240         IF(lwp) WRITE(numout,*) '      quadratic bottom friction'
241         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_bfri2
242         IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_bfri2_max  = ', rn_bfri2_max
243         IF(lwp) WRITE(numout,*) '      background tke   rn_bfeb2  = ', rn_bfeb2
244         IF(lwp) WRITE(numout,*) '      log formulation   ln_bfr2d = ', ln_loglayer
245         IF(lwp) WRITE(numout,*) '      bottom roughness  rn_bfrz0 [m] = ', rn_bfrz0
246         IF( rn_bfrz0<=0.e0 ) THEN
247            WRITE(ctmp1,*) '      bottom roughness must be strictly positive'
248            CALL ctl_stop( ctmp1 )
249         ENDIF
250         IF( ln_bfr2d ) THEN
251            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d
252            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
253         ENDIF
254         !
255         IF(ln_bfr2d) THEN
256            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
257            CALL iom_open('bfr_coef.nc',inum)
258            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
259            CALL iom_close(inum)
260            !
261            bfrcoef2d(:,:) = rn_bfri2 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
262         ELSE
263            bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
264         ENDIF
265         !
266         IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all
267#  if defined key_vectopt_loop
268            DO jj = 1, 1
269!CDIR NOVERRCHK
270               DO ji = 1, jpij   ! vector opt. (forced unrolling)
271#  else
272!CDIR NOVERRCHK
273            DO jj = 1, jpj
274!CDIR NOVERRCHK
275               DO ji = 1, jpi
276#  endif
277                  ikbt = mbkt(ji,jj)
278                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
279                  bfrcoef2d(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
280                  bfrcoef2d(ji,jj) = MIN(bfrcoef2d(ji,jj), rn_bfri2_max)
281               END DO
282            END DO
283         ENDIF
284         !
285      CASE DEFAULT
286         IF(lwp) WRITE(ctmp1,*) '         bad flag value for nn_bfr = ', nn_bfr
287         CALL ctl_stop( ctmp1 )
288         !
289      END SELECT
290      !
291      IF(lwp) WRITE(numout,*) '      implicit bottom friction switch                ln_bfrimp  = ', ln_bfrimp
292      !
293      !                              ! Make sure ln_zdfexp=.false. when use implicit bfr
294      IF( ln_bfrimp .AND. ln_zdfexp ) THEN
295         IF(lwp) THEN
296            WRITE(numout,*)
297            WRITE(numout,*) 'Implicit bottom friction can only be used when ln_zdfexp=.false.'
298            WRITE(numout,*) '         but you set: ln_bfrimp=.true. and ln_zdfexp=.true.!!!!'
299            WRITE(ctmp1,*)  '         set either ln_zdfexp = .false or ln_bfrimp = .false.'
300            CALL ctl_stop( ctmp1 )
301         END IF
302      END IF
303      !
304      ! Basic stability check on bottom friction coefficient
305      !
306      ictu = 0               ! counter for stability criterion breaches at U-pts
307      ictv = 0               ! counter for stability criterion breaches at V-pts
308      zminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient
309      zmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
310      !
311#  if defined key_vectopt_loop
312      DO jj = 1, 1
313!CDIR NOVERRCHK
314         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
315#  else
316!CDIR NOVERRCHK
317      DO jj = 2, jpjm1
318!CDIR NOVERRCHK
319         DO ji = 2, jpim1
320#  endif
321             ikbu = mbku(ji,jj)       ! deepest ocean level at u- and v-points
322             ikbv = mbkv(ji,jj)
323             zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt
324             zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt
325             IF( ABS( bfrcoef2d(ji,jj) ) > zfru ) THEN
326                IF( ln_ctl ) THEN
327                   WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu
328                   WRITE(numout,*) 'BFR ', ABS( bfrcoef2d(ji,jj) ), zfru
329                ENDIF
330                ictu = ictu + 1
331             ENDIF
332             IF( ABS( bfrcoef2d(ji,jj) ) > zfrv ) THEN
333                 IF( ln_ctl ) THEN
334                     WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv
335                     WRITE(numout,*) 'BFR ', bfrcoef2d(ji,jj), zfrv
336                 ENDIF
337                 ictv = ictv + 1
338             ENDIF
339             zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  )
340             zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  )
341         END DO
342      END DO
343      IF( lk_mpp ) THEN
344         CALL mpp_sum( ictu )
345         CALL mpp_sum( ictv )
346         CALL mpp_min( zminbfr )
347         CALL mpp_max( zmaxbfr )
348      ENDIF
349      IF( .NOT.ln_bfrimp) THEN
350      IF( lwp .AND. ictu + ictv > 0 ) THEN
351         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
352         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
353         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
354         WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
355      ENDIF
356      ENDIF
357      !
358      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init')
359      !
360   END SUBROUTINE zdf_bfr_init
361
362   !!======================================================================
363END MODULE zdfbfr
Note: See TracBrowser for help on using the repository browser.