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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyndmp.F90 @ 11001

Last change on this file since 11001 was 11001, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : C1D and step.F90. Passes SETTE. Compiles with key_c1d.
  • Property svn:keywords set to Id
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   !
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      REWIND( numnam_ref )              ! Namelist namc1d_dyndmp in reference namelist :
82      READ  ( numnam_ref, namc1d_dyndmp, IOSTAT = ios, ERR = 901)
83901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namc1d_dyndmp in reference namelist', lwp )
84      REWIND( numnam_cfg )              ! Namelist namc1d_dyndmp in configuration namelist : Parameters of the run
85      READ  ( numnam_cfg, namc1d_dyndmp, IOSTAT = ios, ERR = 902 )
86902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namc1d_dyndmp in configuration namelist', lwp )
87      IF(lwm) WRITE ( numond, namc1d_dyndmp )
88      !
89      IF(lwp) THEN                           ! control print
90         WRITE(numout,*)
91         WRITE(numout,*) 'dyn_dmp_init : U and V current Newtonian damping'
92         WRITE(numout,*) '~~~~~~~~~~~~'
93         WRITE(numout,*) '   Namelist namc1d_dyndmp : Set damping flag'
94         WRITE(numout,*) '      add a damping term or not       ln_dyndmp = ', ln_dyndmp
95         WRITE(numout,*) '   Namelist namtra_dmp    : Set damping parameters'
96         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
97         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp
98         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto
99         WRITE(numout,*)
100      ENDIF
101      !
102      IF( ln_dyndmp ) THEN
103         !                                   !==   allocate the data arrays   ==!
104         IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' )
105         !
106         SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==!
107         CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column'
108         CASE ( 1    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the turbocline (avt > 5 cm2/s)'
109         CASE ( 2    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the mixed layer'
110         CASE DEFAULT
111            WRITE(ctmp1,*) '          bad flag value for nn_zdmp = ', nn_zdmp
112            CALL ctl_stop(ctmp1)
113         END SELECT
114         !
115         IF( .NOT. ln_uvd_dyndmp ) THEN      ! force the initialization of U & V current data for damping
116            CALL ctl_warn( 'dyn_dmp_init: U & V current read data not initialized, we force ln_uvd_dyndmp=T' )
117            CALL dta_uvd_init( ld_dyndmp=ln_dyndmp )
118         ENDIF
119         !
120         utrdmp(:,:,:) = 0._wp               ! internal damping trends
121         vtrdmp(:,:,:) = 0._wp
122         !
123         !Read in mask from file
124         CALL iom_open ( cn_resto, imask)
125         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
126         CALL iom_close( imask )
127      ENDIF
128      !
129   END SUBROUTINE dyn_dmp_init
130
131
132   SUBROUTINE dyn_dmp( kt, Kbb, Kmm, puu, pvv, Krhs )
133      !!----------------------------------------------------------------------
134      !!                   ***  ROUTINE dyn_dmp  ***
135      !!                 
136      !! ** Purpose :   Compute the momentum trends due to a newtonian damping
137      !!      of the ocean velocities towards the given data and add it to the
138      !!      general momentum trends.
139      !!
140      !! ** Method  :   Compute Newtonian damping towards u_dta and v_dta
141      !!      and add to the general momentum trends:
142      !!                     puu(Krhs) = puu(Krhs) + resto_uv * (u_dta - puu(Kbb))
143      !!                     pvv(Krhs) = pvv(Krhs) + resto_uv * (v_dta - pvv(Kbb))
144      !!      The trend is computed either throughout the water column
145      !!      (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or
146      !!      below the well mixed layer (nn_zdmp=2)
147      !!
148      !! ** Action  : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs))   momentum trends updated with the damping trend
149      !!----------------------------------------------------------------------
150      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time-step index
151      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm, Krhs ! ocean time level indices
152      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv       ! ocean velocities and RHS of momentum equation
153      !!
154      INTEGER  ::   ji, jj, jk   ! dummy loop indices
155      REAL(wp) ::   zua, zva     ! local scalars
156      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zuv_dta   ! Read in data
157      !!----------------------------------------------------------------------
158      !
159      IF( ln_timing )   CALL timing_start( 'dyn_dmp' )
160      !
161      !
162      !                           !==   read and interpolate U & V current data at kt   ==!
163      CALL dta_uvd( kt, Kmm, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside
164                                  !!!       the C1D context (use of U,V grid variables)
165      !
166      SELECT CASE ( nn_zdmp )     !==   Calculate/add Newtonian damping to the momentum trend   ==!
167      !
168      CASE( 0 )                   ! Newtonian damping throughout the water column
169         DO jk = 1, jpkm1
170            DO jj = 2, jpjm1
171               DO ji = fs_2, fs_jpim1   ! vector opt.
172                  zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
173                  zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
174                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
175                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
176                  utrdmp(ji,jj,jk) = zua           ! save the trends
177                  vtrdmp(ji,jj,jk) = zva     
178               END DO
179            END DO
180         END DO
181         !
182      CASE ( 1 )                  ! no damping above the turbocline (avt > 5 cm2/s)
183         DO jk = 1, jpkm1
184            DO jj = 2, jpjm1
185               DO ji = fs_2, fs_jpim1   ! vector opt.
186                  IF( avt(ji,jj,jk) <= avt_c ) THEN
187                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
188                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
189                  ELSE
190                     zua = 0._wp
191                     zva = 0._wp 
192                  ENDIF
193                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
194                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
195                  utrdmp(ji,jj,jk) = zua           ! save the trends
196                  vtrdmp(ji,jj,jk) = zva
197               END DO
198            END DO
199         END DO
200         !
201      CASE ( 2 )                  ! no damping in the mixed layer
202         DO jk = 1, jpkm1
203            DO jj = 2, jpjm1
204               DO ji = fs_2, fs_jpim1   ! vector opt.
205                  IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
206                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
207                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
208                  ELSE
209                     zua = 0._wp
210                     zva = 0._wp 
211                  ENDIF
212                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
213                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
214                  utrdmp(ji,jj,jk) = zua           ! save the trends
215                  vtrdmp(ji,jj,jk) = zva
216               END DO
217            END DO
218         END DO
219         !
220      END SELECT
221      !
222      !                           ! Control print
223      IF( ln_ctl   )   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' dmp  - Ua: ', mask1=umask,   &
224         &                           tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
225      !
226      !
227      IF( ln_timing )   CALL timing_stop( 'dyn_dmp')
228      !
229   END SUBROUTINE dyn_dmp
230
231   !!======================================================================
232END MODULE dyndmp
Note: See TracBrowser for help on using the repository browser.