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.
dyndmp.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/C1D – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 11.3 KB
Line 
1MODULE dyndmp
2   !!======================================================================
3   !!                       ***  MODULE  dyndmp  ***
4   !! Ocean dynamics: internal restoring trend on momentum (U and V current)
5   !!                 This should only be used for C1D case in current form
6   !!======================================================================
7   !! History :  3.5  ! 2013-08  (D. Calvert)  Original code
8   !!            3.6  ! 2014-08  (T. Graham) Modified to use netcdf file of
9   !!                             restoration coefficients supplied to tradmp
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   dyn_dmp_alloc : allocate dyndmp arrays
14   !!   dyn_dmp_init  : namelist read, parameter control and resto coeff.
15   !!   dyn_dmp       : update the momentum trend with the internal damping
16   !!----------------------------------------------------------------------
17   USE oce            ! ocean: variables
18   USE dom_oce        ! ocean: domain variables
19   USE c1d            ! 1D vertical configuration
20   USE tradmp         ! ocean: internal damping
21   USE zdf_oce        ! ocean: vertical physics
22   USE phycst         ! physical constants
23   USE dtauvd         ! data: U & V current
24   USE zdfmxl         ! vertical physics: mixed layer depth
25   USE in_out_manager ! I/O manager
26   USE lib_mpp        ! MPP library
27   USE prtctl         ! Print control
28   USE wrk_nemo       ! Memory allocation
29   USE timing         ! Timing
30   USE iom            ! I/O manager
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_dmp_init ! routine called by nemogcm.F90
36   PUBLIC   dyn_dmp      ! routine called by step_c1d.F90
37
38   LOGICAL, PUBLIC ::   ln_dyndmp           ! Flag for Newtonian damping
39
40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  utrdmp    ! damping U current trend (m/s2)
41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  vtrdmp    ! damping V current trend (m/s2)
42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  resto_uv  ! restoring coeff. on U & V current
43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   INTEGER FUNCTION dyn_dmp_alloc()
55      !!----------------------------------------------------------------------
56      !!                ***  FUNCTION dyn_dmp_alloc  ***
57      !!----------------------------------------------------------------------
58      ALLOCATE( utrdmp(jpi,jpj,jpk), vtrdmp(jpi,jpj,jpk), resto_uv(jpi,jpj,jpk), STAT= dyn_dmp_alloc )
59      !
60      IF( lk_mpp            )   CALL mpp_sum ( dyn_dmp_alloc )
61      IF( dyn_dmp_alloc > 0 )   CALL ctl_warn('dyn_dmp_alloc: allocation of arrays failed')
62      !
63   END FUNCTION dyn_dmp_alloc
64
65
66   SUBROUTINE dyn_dmp_init
67      !!----------------------------------------------------------------------
68      !!                  ***  ROUTINE dyn_dmp_init  ***
69      !!
70      !! ** Purpose :   Initialization for the Newtonian damping
71      !!
72      !! ** Method  : - read the ln_dyndmp parameter from the namc1d_dyndmp namelist
73      !!              - allocate damping arrays
74      !!              - check the parameters of the namtra_dmp namelist
75      !!              - calculate damping coefficient
76      !!----------------------------------------------------------------------
77      NAMELIST/namc1d_dyndmp/ ln_dyndmp
78      INTEGER :: ios
79      INTEGER :: imask
80      !!----------------------------------------------------------------------
81
82      REWIND( numnam_ref )              ! Namelist namc1d_dyndmp in reference namelist :
83      READ  ( numnam_ref, namc1d_dyndmp, IOSTAT = ios, ERR = 901)
84901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in reference namelist', lwp )
85
86      REWIND( numnam_cfg )              ! Namelist namc1d_dyndmp in configuration namelist : Parameters of the run
87      READ  ( numnam_cfg, namc1d_dyndmp, IOSTAT = ios, ERR = 902 )
88902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in configuration namelist', lwp )
89      IF(lwm .AND. nprint > 2) WRITE ( numond, namc1d_dyndmp )
90
91      IF(lwp) THEN                           ! control print
92         WRITE(numout,*)
93         WRITE(numout,*) 'dyn_dmp_init : U and V current Newtonian damping'
94         WRITE(numout,*) '~~~~~~~~~~~~'
95         WRITE(numout,*) '   Namelist namc1d_dyndmp : Set damping flag'
96         WRITE(numout,*) '      add a damping term or not       ln_dyndmp = ', ln_dyndmp
97         WRITE(numout,*) '   Namelist namtra_dmp    : Set damping parameters'
98         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
99         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp
100         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto
101         WRITE(numout,*)
102         IF(lflush) CALL flush(numout)
103      ENDIF
104
105      IF( ln_dyndmp ) THEN
106         !                                   !==   allocate the data arrays   ==!
107         IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' )
108         !
109         SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==!
110         CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column'
111         CASE ( 1    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the turbocline (avt > 5 cm2/s)'
112         CASE ( 2    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the mixed layer'
113         CASE DEFAULT
114            WRITE(ctmp1,*) '          bad flag value for nn_zdmp = ', nn_zdmp
115            CALL ctl_stop(ctmp1)
116         END SELECT
117         !
118         IF( .NOT. ln_uvd_dyndmp ) THEN      ! force the initialization of U & V current data for damping
119            CALL ctl_warn( 'dyn_dmp_init: U & V current read data not initialized, we force ln_uvd_dyndmp=T' )
120            CALL dta_uvd_init( ld_dyndmp=ln_dyndmp )
121         ENDIF
122         !
123         utrdmp(:,:,:) = 0._wp               ! internal damping trends
124         vtrdmp(:,:,:) = 0._wp
125         !
126         !Read in mask from file
127         CALL iom_open ( cn_resto, imask)
128         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
129         CALL iom_close( imask )
130      ENDIF
131      !
132      IF(lwp .AND. lflush) CALL flush(numout)
133      !
134   END SUBROUTINE dyn_dmp_init
135
136
137   SUBROUTINE dyn_dmp( kt )
138      !!----------------------------------------------------------------------
139      !!                   ***  ROUTINE dyn_dmp  ***
140      !!                 
141      !! ** Purpose :   Compute the momentum trends due to a newtonian damping
142      !!      of the ocean velocities towards the given data and add it to the
143      !!      general momentum trends.
144      !!
145      !! ** Method  :   Compute Newtonian damping towards u_dta and v_dta
146      !!      and add to the general momentum trends:
147      !!                     ua = ua + resto_uv * (u_dta - ub)
148      !!                     va = va + resto_uv * (v_dta - vb)
149      !!      The trend is computed either throughout the water column
150      !!      (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or
151      !!      below the well mixed layer (nn_zdmp=2)
152      !!
153      !! ** Action  : - (ua,va)   momentum trends updated with the damping trend
154      !!----------------------------------------------------------------------
155      INTEGER, INTENT(in) ::   kt                        ! ocean time-step index
156      !!
157      INTEGER  ::   ji, jj, jk                           ! dummy loop indices
158      REAL(wp) ::   zua, zva                             ! local scalars
159      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zuv_dta  ! Read in data
160      !!----------------------------------------------------------------------
161      !
162      IF( nn_timing == 1 )  CALL timing_start( 'dyn_dmp' )
163      !
164      CALL wrk_alloc( jpi, jpj, jpk, 2,  zuv_dta )
165      !
166      !                           !==   read and interpolate U & V current data at kt   ==!
167      CALL dta_uvd( kt, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside
168                                  !!!       the C1D context (use of U,V grid variables)
169      !
170      SELECT CASE ( nn_zdmp )     !==   Calculate/add Newtonian damping to the momentum trend   ==!
171      !
172      CASE( 0 )                   ! Newtonian damping throughout the water column
173         DO jk = 1, jpkm1
174            DO jj = 2, jpjm1
175               DO ji = fs_2, fs_jpim1   ! vector opt.
176                  zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
177                  zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
178                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
179                  va(ji,jj,jk) = va(ji,jj,jk) + zva
180                  utrdmp(ji,jj,jk) = zua           ! save the trends
181                  vtrdmp(ji,jj,jk) = zva     
182               END DO
183            END DO
184         END DO
185         !
186      CASE ( 1 )                  ! no damping above the turbocline (avt > 5 cm2/s)
187         DO jk = 1, jpkm1
188            DO jj = 2, jpjm1
189               DO ji = fs_2, fs_jpim1   ! vector opt.
190                  IF( avt(ji,jj,jk) <= 5.e-4_wp ) THEN
191                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
192                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
193                  ELSE
194                     zua = 0._wp
195                     zva = 0._wp 
196                  ENDIF
197                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
198                  va(ji,jj,jk) = va(ji,jj,jk) + zva
199                  utrdmp(ji,jj,jk) = zua           ! save the trends
200                  vtrdmp(ji,jj,jk) = zva
201               END DO
202            END DO
203         END DO
204         !
205      CASE ( 2 )                  ! no damping in the mixed layer
206         DO jk = 1, jpkm1
207            DO jj = 2, jpjm1
208               DO ji = fs_2, fs_jpim1   ! vector opt.
209                  IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
210                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
211                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
212                  ELSE
213                     zua = 0._wp
214                     zva = 0._wp 
215                  ENDIF
216                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
217                  va(ji,jj,jk) = va(ji,jj,jk) + zva
218                  utrdmp(ji,jj,jk) = zua           ! save the trends
219                  vtrdmp(ji,jj,jk) = zva
220               END DO
221            END DO
222         END DO
223         !
224      END SELECT
225      !
226!!gm      !                           ! Trend diagnostic
227!!gm      IF( l_trddyn )   CALL trd_mod( utrdmp, vtrdmp, jpdyn_trd_dat, 'DYN', kt )
228      !
229      !                           ! Control print
230      IF( ln_ctl   )   CALL prt_ctl( tab3d_1=ua(:,:,:), clinfo1=' dmp  - Ua: ', mask1=umask,   &
231         &                           tab3d_2=va(:,:,:), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
232      !
233      CALL wrk_dealloc( jpi, jpj, jpk, 2,  zuv_dta )
234      !
235      IF( nn_timing == 1 )  CALL timing_stop( 'dyn_dmp')
236      !
237   END SUBROUTINE dyn_dmp
238
239   !!======================================================================
240END MODULE dyndmp
Note: See TracBrowser for help on using the repository browser.