-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdyndmp.F90
237 lines (221 loc) · 11.2 KB
/
dyndmp.F90
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
151
152
153
154
155
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
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
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 wrk_nemo ! Memory allocation
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 "domzgr_substitute.h90"
# include "vectopt_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OPA 3.3 , NEMO Consortium (2010)
!! $Id$
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
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 )
!
IF( lk_mpp ) CALL mpp_sum ( 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
!!----------------------------------------------------------------------
NAMELIST/namc1d_dyndmp/ ln_dyndmp
INTEGER :: ios
INTEGER :: imask
!!----------------------------------------------------------------------
REWIND( numnam_ref ) ! Namelist namc1d_dyndmp in reference namelist :
READ ( numnam_ref, namc1d_dyndmp, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in reference namelist', lwp )
REWIND( numnam_cfg ) ! Namelist namc1d_dyndmp in configuration namelist : Parameters of the run
READ ( numnam_cfg, namc1d_dyndmp, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in configuration namelist', lwp )
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_autoglo, 'resto', resto)
CALL iom_close( imask )
ENDIF
!
END SUBROUTINE dyn_dmp_init
SUBROUTINE dyn_dmp( kt )
!!----------------------------------------------------------------------
!! *** 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:
!! ua = ua + resto_uv * (u_dta - ub)
!! va = va + resto_uv * (v_dta - vb)
!! 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 : - (ua,va) momentum trends updated with the damping trend
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step index
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zua, zva ! local scalars
REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zuv_dta ! Read in data
!!----------------------------------------------------------------------
!
IF( nn_timing == 1 ) CALL timing_start( 'dyn_dmp' )
!
CALL wrk_alloc( jpi, jpj, jpk, 2, zuv_dta )
!
! !== read and interpolate U & V current data at kt ==!
CALL dta_uvd( kt, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside
!!! the C1D context (use of U,V grid variables)
!
SELECT CASE ( nn_zdmp ) !== Calculate/add Newtonian damping to the momentum trend ==!
!
CASE( 0 ) ! Newtonian damping throughout the water column
DO jk = 1, jpkm1
DO jj = 2, jpjm1
DO ji = fs_2, fs_jpim1 ! vector opt.
zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
ua(ji,jj,jk) = ua(ji,jj,jk) + zua
va(ji,jj,jk) = va(ji,jj,jk) + zva
utrdmp(ji,jj,jk) = zua ! save the trends
vtrdmp(ji,jj,jk) = zva
END DO
END DO
END DO
!
CASE ( 1 ) ! no damping above the turbocline (avt > 5 cm2/s)
DO jk = 1, jpkm1
DO jj = 2, jpjm1
DO ji = fs_2, fs_jpim1 ! vector opt.
IF( avt(ji,jj,jk) <= 5.e-4_wp ) THEN
zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
ELSE
zua = 0._wp
zva = 0._wp
ENDIF
ua(ji,jj,jk) = ua(ji,jj,jk) + zua
va(ji,jj,jk) = va(ji,jj,jk) + zva
utrdmp(ji,jj,jk) = zua ! save the trends
vtrdmp(ji,jj,jk) = zva
END DO
END DO
END DO
!
CASE ( 2 ) ! no damping in the mixed layer
DO jk = 1, jpkm1
DO jj = 2, jpjm1
DO ji = fs_2, fs_jpim1 ! vector opt.
IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
ELSE
zua = 0._wp
zva = 0._wp
ENDIF
ua(ji,jj,jk) = ua(ji,jj,jk) + zua
va(ji,jj,jk) = va(ji,jj,jk) + zva
utrdmp(ji,jj,jk) = zua ! save the trends
vtrdmp(ji,jj,jk) = zva
END DO
END DO
END DO
!
END SELECT
!
!!gm ! ! Trend diagnostic
!!gm IF( l_trddyn ) CALL trd_mod( utrdmp, vtrdmp, jpdyn_trd_dat, 'DYN', kt )
!
! ! Control print
IF( ln_ctl ) CALL prt_ctl( tab3d_1=ua(:,:,:), clinfo1=' dmp - Ua: ', mask1=umask, &
& tab3d_2=va(:,:,:), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
CALL wrk_dealloc( jpi, jpj, jpk, 2, zuv_dta )
!
IF( nn_timing == 1 ) CALL timing_stop( 'dyn_dmp')
!
END SUBROUTINE dyn_dmp
!!======================================================================
END MODULE dyndmp