Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
MODULE dyndmp
!!======================================================================
!! *** MODULE dyndmp ***
!! Ocean dynamics: internal restoring trend on momentum (U and V current)
!! This should only be used for C1D case in current form
!!======================================================================
!! History : 3.5 ! 2013-08 (D. Calvert) Original code
!! 3.6 ! 2014-08 (T. Graham) Modified to use netcdf file of
!! restoration coefficients supplied to tradmp
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! dyn_dmp_alloc : allocate dyndmp arrays
!! dyn_dmp_init : namelist read, parameter control and resto coeff.
!! dyn_dmp : update the momentum trend with the internal damping
!!----------------------------------------------------------------------
USE oce ! ocean: variables
USE dom_oce ! ocean: domain variables
USE c1d ! 1D vertical configuration
USE tradmp ! ocean: internal damping
USE zdf_oce ! ocean: vertical physics
USE phycst ! physical constants
USE dtauvd ! data: U & V current
USE zdfmxl ! vertical physics: mixed layer depth
!
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
USE prtctl ! Print control
USE timing ! Timing
USE iom ! I/O manager
IMPLICIT NONE
PRIVATE
PUBLIC dyn_dmp_init ! routine called by nemogcm.F90
PUBLIC dyn_dmp ! routine called by step_c1d.F90
LOGICAL, PUBLIC :: ln_dyndmp !: Flag for Newtonian damping
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: utrdmp !: damping U current trend (m/s2)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vtrdmp !: damping V current trend (m/s2)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: resto_uv !: restoring coeff. on U & V current
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dyndmp.F90 15062 2021-06-28 11:19:48Z jchanut $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
INTEGER FUNCTION dyn_dmp_alloc()
!!----------------------------------------------------------------------
!! *** FUNCTION dyn_dmp_alloc ***
!!----------------------------------------------------------------------
ALLOCATE( utrdmp(jpi,jpj,jpk), vtrdmp(jpi,jpj,jpk), resto_uv(jpi,jpj,jpk), STAT= dyn_dmp_alloc )
!
CALL mpp_sum ( 'dyndmp', dyn_dmp_alloc )
IF( dyn_dmp_alloc > 0 ) CALL ctl_warn('dyn_dmp_alloc: allocation of arrays failed')
!
END FUNCTION dyn_dmp_alloc
SUBROUTINE dyn_dmp_init
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_dmp_init ***
!!
!! ** Purpose : Initialization for the Newtonian damping
!!
!! ** Method : - read the ln_dyndmp parameter from the namc1d_dyndmp namelist
!! - allocate damping arrays
!! - check the parameters of the namtra_dmp namelist
!! - calculate damping coefficient
!!----------------------------------------------------------------------
INTEGER :: ios, imask ! local integers
!!
NAMELIST/namc1d_dyndmp/ ln_dyndmp
!!----------------------------------------------------------------------
!
READ ( numnam_ref, namc1d_dyndmp, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in reference namelist' )
READ ( numnam_cfg, namc1d_dyndmp, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in configuration namelist' )
IF(lwm) WRITE ( numond, namc1d_dyndmp )
!
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) 'dyn_dmp_init : U and V current Newtonian damping'
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namc1d_dyndmp : Set damping flag'
WRITE(numout,*) ' add a damping term or not ln_dyndmp = ', ln_dyndmp
WRITE(numout,*) ' Namelist namtra_dmp : Set damping parameters'
WRITE(numout,*) ' Apply relaxation or not ln_tradmp = ', ln_tradmp
WRITE(numout,*) ' mixed layer damping option nn_zdmp = ', nn_zdmp
WRITE(numout,*) ' Damping file name cn_resto = ', cn_resto
WRITE(numout,*)
ENDIF
!
IF( ln_dyndmp ) THEN
! !== allocate the data arrays ==!
IF( dyn_dmp_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' )
!
SELECT CASE ( nn_zdmp ) !== control print of vertical option ==!
CASE ( 0 ) ; IF(lwp) WRITE(numout,*) ' momentum damping throughout the water column'
CASE ( 1 ) ; IF(lwp) WRITE(numout,*) ' no momentum damping in the turbocline (avt > 5 cm2/s)'
CASE ( 2 ) ; IF(lwp) WRITE(numout,*) ' no momentum damping in the mixed layer'
CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for nn_zdmp = ', nn_zdmp
CALL ctl_stop(ctmp1)
END SELECT
!
IF( .NOT. ln_uvd_dyndmp ) THEN ! force the initialization of U & V current data for damping
CALL ctl_warn( 'dyn_dmp_init: U & V current read data not initialized, we force ln_uvd_dyndmp=T' )
CALL dta_uvd_init( ld_dyndmp=ln_dyndmp )
ENDIF
!
utrdmp(:,:,:) = 0._wp ! internal damping trends
vtrdmp(:,:,:) = 0._wp
!
!Read in mask from file
CALL iom_open ( cn_resto, imask)
CALL iom_get ( imask, jpdom_auto, 'resto', resto)
CALL iom_close( imask )
ENDIF
!
END SUBROUTINE dyn_dmp_init
SUBROUTINE dyn_dmp( kt, Kbb, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_dmp ***
!!
!! ** Purpose : Compute the momentum trends due to a newtonian damping
!! of the ocean velocities towards the given data and add it to the
!! general momentum trends.
!!
!! ** Method : Compute Newtonian damping towards u_dta and v_dta
!! and add to the general momentum trends:
!! puu(Krhs) = puu(Krhs) + resto_uv * (u_dta - puu(Kbb))
!! pvv(Krhs) = pvv(Krhs) + resto_uv * (v_dta - pvv(Kbb))
!! The trend is computed either throughout the water column
!! (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or
!! below the well mixed layer (nn_zdmp=2)
!!
!! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) momentum trends updated with the damping trend
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices
sparonuz
committed
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zua, zva ! local scalars
sparonuz
committed
REAL(dp), DIMENSION(jpi,jpj,jpk,2) :: zuv_dta ! Read in data
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start( 'dyn_dmp' )
!
!
! !== read and interpolate U & V current data at kt ==!
CALL dta_uvd( kt, Kmm, zuv_dta(:,:,:,1), zuv_dta(:,:,:,2))
!
SELECT CASE ( nn_zdmp ) !== Calculate/add Newtonian damping to the momentum trend ==!
!
CASE( 0 ) ! Newtonian damping throughout the water column
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
utrdmp(ji,jj,jk) = zua ! save the trends
vtrdmp(ji,jj,jk) = zva
END_3D
!
CASE ( 1 ) ! no damping above the turbocline (avt > 5 cm2/s)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
IF( avt(ji,jj,jk) <= avt_c ) THEN
zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
ELSE
zua = 0._wp
zva = 0._wp
ENDIF
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
utrdmp(ji,jj,jk) = zua ! save the trends
vtrdmp(ji,jj,jk) = zva
END_3D
!
CASE ( 2 ) ! no damping in the mixed layer
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
ELSE
zua = 0._wp
zva = 0._wp
ENDIF
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
utrdmp(ji,jj,jk) = zua ! save the trends
vtrdmp(ji,jj,jk) = zva
END_3D
!
END SELECT
!
! ! Control print
IF( sn_cfctl%l_prtctl ) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' dmp - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
!
IF( ln_timing ) CALL timing_stop( 'dyn_dmp')
!
END SUBROUTINE dyn_dmp
!!======================================================================
END MODULE dyndmp