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 NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/C1D – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/C1D/dyndmp.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 10.6 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   !
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! MPP library
28   USE prtctl         ! Print control
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 "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
48   !! $Id$
49   !! Software governed by the CeCILL license (see ./LICENSE)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   INTEGER FUNCTION dyn_dmp_alloc()
54      !!----------------------------------------------------------------------
55      !!                ***  FUNCTION dyn_dmp_alloc  ***
56      !!----------------------------------------------------------------------
57      ALLOCATE( utrdmp(jpi,jpj,jpk), vtrdmp(jpi,jpj,jpk), resto_uv(jpi,jpj,jpk), STAT= dyn_dmp_alloc )
58      !
59      CALL mpp_sum ( 'dyndmp', dyn_dmp_alloc )
60      IF( dyn_dmp_alloc > 0 )   CALL ctl_warn('dyn_dmp_alloc: allocation of arrays failed')
61      !
62   END FUNCTION dyn_dmp_alloc
63
64
65   SUBROUTINE dyn_dmp_init
66      !!----------------------------------------------------------------------
67      !!                  ***  ROUTINE dyn_dmp_init  ***
68      !!
69      !! ** Purpose :   Initialization for the Newtonian damping
70      !!
71      !! ** Method  : - read the ln_dyndmp parameter from the namc1d_dyndmp namelist
72      !!              - allocate damping arrays
73      !!              - check the parameters of the namtra_dmp namelist
74      !!              - calculate damping coefficient
75      !!----------------------------------------------------------------------
76      INTEGER ::   ios, imask   ! local integers
77      !!
78      NAMELIST/namc1d_dyndmp/ ln_dyndmp
79      !!----------------------------------------------------------------------
80      !
81      READ  ( numnam_ref, namc1d_dyndmp, IOSTAT = ios, ERR = 901)
82901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namc1d_dyndmp in reference namelist' )
83      READ  ( numnam_cfg, namc1d_dyndmp, IOSTAT = ios, ERR = 902 )
84902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namc1d_dyndmp in configuration namelist' )
85      IF(lwm) WRITE ( numond, namc1d_dyndmp )
86      !
87      IF(lwp) THEN                           ! control print
88         WRITE(numout,*)
89         WRITE(numout,*) 'dyn_dmp_init : U and V current Newtonian damping'
90         WRITE(numout,*) '~~~~~~~~~~~~'
91         WRITE(numout,*) '   Namelist namc1d_dyndmp : Set damping flag'
92         WRITE(numout,*) '      add a damping term or not       ln_dyndmp = ', ln_dyndmp
93         WRITE(numout,*) '   Namelist namtra_dmp    : Set damping parameters'
94         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
95         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp
96         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto
97         WRITE(numout,*)
98      ENDIF
99      !
100      IF( ln_dyndmp ) THEN
101         !                                   !==   allocate the data arrays   ==!
102         IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' )
103         !
104         SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==!
105         CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column'
106         CASE ( 1    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the turbocline (avt > 5 cm2/s)'
107         CASE ( 2    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the mixed layer'
108         CASE DEFAULT
109            WRITE(ctmp1,*) '          bad flag value for nn_zdmp = ', nn_zdmp
110            CALL ctl_stop(ctmp1)
111         END SELECT
112         !
113         IF( .NOT. ln_uvd_dyndmp ) THEN      ! force the initialization of U & V current data for damping
114            CALL ctl_warn( 'dyn_dmp_init: U & V current read data not initialized, we force ln_uvd_dyndmp=T' )
115            CALL dta_uvd_init( ld_dyndmp=ln_dyndmp )
116         ENDIF
117         !
118         utrdmp(:,:,:) = 0._wp               ! internal damping trends
119         vtrdmp(:,:,:) = 0._wp
120         !
121         !Read in mask from file
122         CALL iom_open ( cn_resto, imask)
123         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
124         CALL iom_close( imask )
125      ENDIF
126      !
127   END SUBROUTINE dyn_dmp_init
128
129
130   SUBROUTINE dyn_dmp( kt )
131      !!----------------------------------------------------------------------
132      !!                   ***  ROUTINE dyn_dmp  ***
133      !!                 
134      !! ** Purpose :   Compute the momentum trends due to a newtonian damping
135      !!      of the ocean velocities towards the given data and add it to the
136      !!      general momentum trends.
137      !!
138      !! ** Method  :   Compute Newtonian damping towards u_dta and v_dta
139      !!      and add to the general momentum trends:
140      !!                     ua = ua + resto_uv * (u_dta - ub)
141      !!                     va = va + resto_uv * (v_dta - vb)
142      !!      The trend is computed either throughout the water column
143      !!      (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or
144      !!      below the well mixed layer (nn_zdmp=2)
145      !!
146      !! ** Action  : - (ua,va)   momentum trends updated with the damping trend
147      !!----------------------------------------------------------------------
148      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
149      !!
150      INTEGER  ::   ji, jj, jk   ! dummy loop indices
151      REAL(wp) ::   zua, zva     ! local scalars
152      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zuv_dta   ! Read in data
153      !!----------------------------------------------------------------------
154      !
155      IF( ln_timing )   CALL timing_start( 'dyn_dmp' )
156      !
157      !
158      !                           !==   read and interpolate U & V current data at kt   ==!
159      CALL dta_uvd( kt, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside
160                                  !!!       the C1D context (use of U,V grid variables)
161      !
162      SELECT CASE ( nn_zdmp )     !==   Calculate/add Newtonian damping to the momentum trend   ==!
163      !
164      CASE( 0 )                   ! Newtonian damping throughout the water column
165         DO jk = 1, jpkm1
166            DO jj = 2, jpjm1
167               DO ji = fs_2, fs_jpim1   ! vector opt.
168                  zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
169                  zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
170                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
171                  va(ji,jj,jk) = va(ji,jj,jk) + zva
172                  utrdmp(ji,jj,jk) = zua           ! save the trends
173                  vtrdmp(ji,jj,jk) = zva     
174               END DO
175            END DO
176         END DO
177         !
178      CASE ( 1 )                  ! no damping above the turbocline (avt > 5 cm2/s)
179         DO jk = 1, jpkm1
180            DO jj = 2, jpjm1
181               DO ji = fs_2, fs_jpim1   ! vector opt.
182                  IF( avt(ji,jj,jk) <= avt_c ) THEN
183                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
184                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
185                  ELSE
186                     zua = 0._wp
187                     zva = 0._wp 
188                  ENDIF
189                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
190                  va(ji,jj,jk) = va(ji,jj,jk) + zva
191                  utrdmp(ji,jj,jk) = zua           ! save the trends
192                  vtrdmp(ji,jj,jk) = zva
193               END DO
194            END DO
195         END DO
196         !
197      CASE ( 2 )                  ! no damping in the mixed layer
198         DO jk = 1, jpkm1
199            DO jj = 2, jpjm1
200               DO ji = fs_2, fs_jpim1   ! vector opt.
201                  IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN
202                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
203                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
204                  ELSE
205                     zua = 0._wp
206                     zva = 0._wp 
207                  ENDIF
208                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
209                  va(ji,jj,jk) = va(ji,jj,jk) + zva
210                  utrdmp(ji,jj,jk) = zua           ! save the trends
211                  vtrdmp(ji,jj,jk) = zva
212               END DO
213            END DO
214         END DO
215         !
216      END SELECT
217      !
218      !                           ! Control print
219      IF( ln_ctl   )   CALL prt_ctl( tab3d_1=ua(:,:,:), clinfo1=' dmp  - Ua: ', mask1=umask,   &
220         &                           tab3d_2=va(:,:,:), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
221      !
222      !
223      IF( ln_timing )   CALL timing_stop( 'dyn_dmp')
224      !
225   END SUBROUTINE dyn_dmp
226
227   !!======================================================================
228END MODULE dyndmp
Note: See TracBrowser for help on using the repository browser.