source: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90 @ 4739

Last change on this file since 4739 was 4739, checked in by timgraham, 7 years ago

Updated C1D/dyndmp.F90 and trcdmp.F90 to read restoration coefficient from a file.
Modified namelist_top_ref to match new options
Bug fixes to DMP_TOOLS tool and addition of custom.F90 to allow users to make modifications. Also changed to use working precision (wp) throughout.

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 trdmod         ! ocean: trends
21   USE trdmod_oce     ! ocean: trends variables
22   USE tradmp         ! ocean: internal damping
23   USE zdf_oce        ! ocean: vertical physics
24   USE phycst         ! physical constants
25   USE dtauvd         ! data: U & V current
26   USE zdfmxl         ! vertical physics: mixed layer depth
27   USE in_out_manager ! I/O manager
28   USE lib_mpp        ! MPP library
29   USE prtctl         ! Print control
30   USE wrk_nemo       ! Memory allocation
31   USE timing         ! Timing
32   USE iom            ! I/O manager
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   dyn_dmp_init ! routine called by nemogcm.F90
38   PUBLIC   dyn_dmp      ! routine called by step_c1d.F90
39
40   LOGICAL, PUBLIC ::   ln_dyndmp           ! Flag for Newtonian damping
41
42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  utrdmp    ! damping U current trend (m/s2)
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  vtrdmp    ! damping V current trend (m/s2)
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  resto_uv  ! restoring coeff. on U & V current
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id: dyndmp.F90 3294 2012-01-28 16:44:18Z rblod $
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   INTEGER FUNCTION dyn_dmp_alloc()
57      !!----------------------------------------------------------------------
58      !!                ***  FUNCTION dyn_dmp_alloc  ***
59      !!----------------------------------------------------------------------
60      ALLOCATE( utrdmp(jpi,jpj,jpk), vtrdmp(jpi,jpj,jpk), resto_uv(jpi,jpj,jpk), STAT= dyn_dmp_alloc )
61      !
62      IF( lk_mpp            )   CALL mpp_sum ( dyn_dmp_alloc )
63      IF( dyn_dmp_alloc > 0 )   CALL ctl_warn('dyn_dmp_alloc: allocation of arrays failed')
64      !
65   END FUNCTION dyn_dmp_alloc
66
67
68   SUBROUTINE dyn_dmp_init
69      !!----------------------------------------------------------------------
70      !!                  ***  ROUTINE dyn_dmp_init  ***
71      !!
72      !! ** Purpose :   Initialization for the Newtonian damping
73      !!
74      !! ** Method  : - read the ln_dyndmp parameter from the namc1d_dyndmp namelist
75      !!              - allocate damping arrays
76      !!              - check the parameters of the namtra_dmp namelist
77      !!              - calculate damping coefficient
78      !!----------------------------------------------------------------------
79      NAMELIST/namc1d_dyndmp/ ln_dyndmp
80      INTEGER :: ios
81      INTEGER :: imask
82      !!----------------------------------------------------------------------
83
84      REWIND( numnam_ref )              ! Namelist namc1d_dyndmp in reference namelist :
85      READ  ( numnam_ref, namc1d_dyndmp, IOSTAT = ios, ERR = 901)
86901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in reference namelist', lwp )
87
88      REWIND( numnam_cfg )              ! Namelist namc1d_dyndmp in configuration namelist : Parameters of the run
89      READ  ( numnam_cfg, namc1d_dyndmp, IOSTAT = ios, ERR = 902 )
90902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in configuration namelist', lwp )
91      IF(lwm) WRITE ( numond, namc1d_dyndmp )
92
93      IF(lwp) THEN                           ! control print
94         WRITE(numout,*)
95         WRITE(numout,*) 'dyn_dmp_init : U and V current Newtonian damping'
96         WRITE(numout,*) '~~~~~~~~~~~~'
97         WRITE(numout,*) '   Namelist namc1d_dyndmp : Set damping flag'
98         WRITE(numout,*) '      add a damping term or not       ln_dyndmp = ', ln_dyndmp
99         WRITE(numout,*) '   Namelist namtra_dmp    : Set damping parameters'
100         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
101         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp
102         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto
103         WRITE(numout,*)
104      ENDIF
105
106      IF( ln_dyndmp ) THEN
107         !                                   !==   allocate the data arrays   ==!
108         IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' )
109         !
110         SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==!
111         CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column'
112         CASE ( 1    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the turbocline (avt > 5 cm2/s)'
113         CASE ( 2    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the mixed layer'
114         CASE DEFAULT
115            WRITE(ctmp1,*) '          bad flag value for nn_zdmp = ', nn_zdmp
116            CALL ctl_stop(ctmp1)
117         END SELECT
118         !
119         IF( .NOT. ln_uvd_dyndmp ) THEN      ! force the initialization of U & V current data for damping
120            CALL ctl_warn( 'dyn_dmp_init: U & V current read data not initialized, we force ln_uvd_dyndmp=T' )
121            CALL dta_uvd_init( ld_dyndmp=ln_dyndmp )
122         ENDIF
123         !
124         utrdmp(:,:,:) = 0._wp               ! internal damping trends
125         vtrdmp(:,:,:) = 0._wp
126         !
127         !Read in mask from file
128         CALL iom_open ( cn_resto, imask)
129         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
130         CALL iom_close( imask )
131      ENDIF
132      !
133   END SUBROUTINE dyn_dmp_init
134
135
136   SUBROUTINE dyn_dmp( kt )
137      !!----------------------------------------------------------------------
138      !!                   ***  ROUTINE dyn_dmp  ***
139      !!                 
140      !! ** Purpose :   Compute the momentum trends due to a newtonian damping
141      !!      of the ocean velocities towards the given data and add it to the
142      !!      general momentum trends.
143      !!
144      !! ** Method  :   Compute Newtonian damping towards u_dta and v_dta
145      !!      and add to the general momentum trends:
146      !!                     ua = ua + resto_uv * (u_dta - ub)
147      !!                     va = va + resto_uv * (v_dta - vb)
148      !!      The trend is computed either throughout the water column
149      !!      (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or
150      !!      below the well mixed layer (nn_zdmp=2)
151      !!
152      !! ** Action  : - (ua,va)   momentum trends updated with the damping trend
153      !!----------------------------------------------------------------------
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      !                           ! Trend diagnostic
227      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.