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 branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/C1D – NEMO

source: branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90 @ 4245

Last change on this file since 4245 was 4245, checked in by cetlod, 10 years ago

dev_locean_cmcc_ingv_ukmo_merc : merge in the MERC_UKMO dev branch with trunk rev 4119
File size: 12.3 KB
Line 
1MODULE dyndmp
2   !!======================================================================
3   !!                       ***  MODULE  dyndmp  ***
4   !! Ocean dynamics: internal restoring trend on momentum (U and V current)
5   !!======================================================================
6   !! History :  3.5  ! 2013-08  (D. Calvert)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   dyn_dmp_alloc : allocate dyndmp arrays
11   !!   dyn_dmp_init  : namelist read, parameter control and resto coeff.
12   !!   dyn_dmp       : update the momentum trend with the internal damping
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean: variables
15   USE dom_oce        ! ocean: domain variables
16   USE c1d            ! 1D vertical configuration
17   USE trdmod         ! ocean: trends
18   USE tradmp         ! ocean: internal damping
19   USE zdf_oce        ! ocean: vertical physics
20   USE phycst         ! physical constants
21   USE dtauvd         ! data: U & V current
22   USE zdfmxl         ! vertical physics: mixed layer depth
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE prtctl         ! Print control
26   USE wrk_nemo       ! Memory allocation
27   USE timing         ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   dyn_dmp_init ! routine called by nemogcm.F90
33   PUBLIC   dyn_dmp      ! routine called by step_c1d.F90
34
35   LOGICAL, PUBLIC ::   ln_dyndmp           ! Flag for Newtonian damping
36
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  utrdmp    ! damping U current trend (m/s2)
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  vtrdmp    ! damping V current trend (m/s2)
39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  resto_uv  ! restoring coeff. on U & V current
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
46   !! $Id: dyndmp.F90 3294 2012-01-28 16:44:18Z rblod $
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   INTEGER FUNCTION dyn_dmp_alloc()
52      !!----------------------------------------------------------------------
53      !!                ***  FUNCTION dyn_dmp_alloc  ***
54      !!----------------------------------------------------------------------
55      ALLOCATE( utrdmp(jpi,jpj,jpk), vtrdmp(jpi,jpj,jpk), resto_uv(jpi,jpj,jpk), STAT= dyn_dmp_alloc )
56      !
57      IF( lk_mpp            )   CALL mpp_sum ( dyn_dmp_alloc )
58      IF( dyn_dmp_alloc > 0 )   CALL ctl_warn('dyn_dmp_alloc: allocation of arrays failed')
59      !
60   END FUNCTION dyn_dmp_alloc
61
62
63   SUBROUTINE dyn_dmp_init
64      !!----------------------------------------------------------------------
65      !!                  ***  ROUTINE dyn_dmp_init  ***
66      !!
67      !! ** Purpose :   Initialization for the Newtonian damping
68      !!
69      !! ** Method  : - read the ln_dyndmp parameter from the namc1d_dyndmp namelist
70      !!              - allocate damping arrays
71      !!              - check the parameters of the namtra_dmp namelist
72      !!              - calculate damping coefficient
73      !!----------------------------------------------------------------------
74      NAMELIST/namc1d_dyndmp/ ln_dyndmp
75      INTEGER :: ios
76      !!----------------------------------------------------------------------
77
78      REWIND( numnam_ref )              ! Namelist namc1d_dyndmp in reference namelist :
79      READ  ( numnam_ref, namc1d_dyndmp, IOSTAT = ios, ERR = 901)
80901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in reference namelist', lwp )
81
82      REWIND( numnam_cfg )              ! Namelist namc1d_dyndmp in configuration namelist : Parameters of the run
83      READ  ( numnam_cfg, namc1d_dyndmp, IOSTAT = ios, ERR = 902 )
84902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d_dyndmp in configuration namelist', lwp )
85      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,*) '      horizontal damping option       nn_hdmp   = ', nn_hdmp
95         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)'
96         WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf
97         WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot
98         WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep
99         WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file
100         WRITE(numout,*)
101      ENDIF
102
103      IF( ln_dyndmp ) THEN
104         !                                   !==   allocate the data arrays   ==!
105         IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' )
106         !
107#if ! defined key_c1d
108         SELECT CASE ( nn_hdmp )             !==   control print of horizontal option   ==!
109         CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   momentum damping in the Med & Red seas only'
110         CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   momentum damping poleward of', nn_hdmp, ' degrees'
111         CASE DEFAULT
112            WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp
113            CALL ctl_stop(ctmp1)
114         END SELECT
115         !
116#endif
117         SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==!
118         CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column'
119         CASE ( 1    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the turbocline (avt > 5 cm2/s)'
120         CASE ( 2    )   ;   IF(lwp) WRITE(numout,*) '   no momentum damping in the mixed layer'
121         CASE DEFAULT
122            WRITE(ctmp1,*) '          bad flag value for nn_zdmp = ', nn_zdmp
123            CALL ctl_stop(ctmp1)
124         END SELECT
125         !
126         IF( .NOT. ln_uvd_dyndmp ) THEN      ! force the initialization of U & V current data for damping
127            CALL ctl_warn( 'dyn_dmp_init: U & V current read data not initialized, we force ln_uvd_dyndmp=T' )
128            CALL dta_uvd_init( ld_dyndmp=ln_dyndmp )
129         ENDIF
130         !
131         utrdmp(:,:,:) = 0._wp               ! internal damping trends
132         vtrdmp(:,:,:) = 0._wp
133         !                                   !==   Damping coefficients calculation:                           ==!
134         !                                   !==   use tradmp.F90 subroutines dtacof, dtacof_zoom and cofdis   ==!
135         !                                   !!!   NOTE: these need to be altered for use in this module if
136         !                                   !!!       they are to be used outside the C1D context
137         !                                   !!!       (use of U,V grid variables)
138         IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto_uv )
139         ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'DYN', resto_uv )
140         ENDIF
141         !
142      ENDIF
143      !
144   END SUBROUTINE dyn_dmp_init
145
146
147   SUBROUTINE dyn_dmp( kt )
148      !!----------------------------------------------------------------------
149      !!                   ***  ROUTINE dyn_dmp  ***
150      !!                 
151      !! ** Purpose :   Compute the momentum trends due to a newtonian damping
152      !!      of the ocean velocities towards the given data and add it to the
153      !!      general momentum trends.
154      !!
155      !! ** Method  :   Compute Newtonian damping towards u_dta and v_dta
156      !!      and add to the general momentum trends:
157      !!                     ua = ua + resto_uv * (u_dta - ub)
158      !!                     va = va + resto_uv * (v_dta - vb)
159      !!      The trend is computed either throughout the water column
160      !!      (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or
161      !!      below the well mixed layer (nn_zdmp=2)
162      !!
163      !! ** Action  : - (ua,va)   momentum trends updated with the damping trend
164      !!----------------------------------------------------------------------
165      !
166      INTEGER, INTENT(in) ::   kt                        ! ocean time-step index
167      !!
168      INTEGER  ::   ji, jj, jk                           ! dummy loop indices
169      REAL(wp) ::   zua, zva                             ! local scalars
170      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zuv_dta  ! Read in data
171      !!----------------------------------------------------------------------
172      !
173      IF( nn_timing == 1 )  CALL timing_start( 'dyn_dmp' )
174      !
175      CALL wrk_alloc( jpi, jpj, jpk, 2,  zuv_dta )
176      !
177      !                           !==   read and interpolate U & V current data at kt   ==!
178      CALL dta_uvd( kt, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside
179                                  !!!       the C1D context (use of U,V grid variables)
180      !
181      SELECT CASE ( nn_zdmp )     !==   Calculate/add Newtonian damping to the momentum trend   ==!
182      !
183      CASE( 0 )                   ! Newtonian damping throughout the water column
184         DO jk = 1, jpkm1
185            DO jj = 2, jpjm1
186               DO ji = fs_2, fs_jpim1   ! vector opt.
187                  zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
188                  zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
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 ( 1 )                  ! no damping above the turbocline (avt > 5 cm2/s)
198         DO jk = 1, jpkm1
199            DO jj = 2, jpjm1
200               DO ji = fs_2, fs_jpim1   ! vector opt.
201                  IF( avt(ji,jj,jk) <= 5.e-4_wp ) 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      CASE ( 2 )                  ! no damping in the mixed layer
217         DO jk = 1, jpkm1
218            DO jj = 2, jpjm1
219               DO ji = fs_2, fs_jpim1   ! vector opt.
220                  IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
221                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
222                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
223                  ELSE
224                     zua = 0._wp
225                     zva = 0._wp 
226                  ENDIF
227                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
228                  va(ji,jj,jk) = va(ji,jj,jk) + zva
229                  utrdmp(ji,jj,jk) = zua           ! save the trends
230                  vtrdmp(ji,jj,jk) = zva
231               END DO
232            END DO
233         END DO
234         !
235      END SELECT
236      !
237      !                           ! Trend diagnostic
238      IF( l_trddyn )   CALL trd_mod( utrdmp, vtrdmp, jpdyn_trd_dat, 'DYN', kt )
239      !
240      !                           ! Control print
241      IF( ln_ctl   )   CALL prt_ctl( tab3d_1=ua(:,:,:), clinfo1=' dmp  - Ua: ', mask1=umask,   &
242         &                           tab3d_2=va(:,:,:), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
243      !
244      CALL wrk_dealloc( jpi, jpj, jpk, 2,  zuv_dta )
245      !
246      IF( nn_timing == 1 )  CALL timing_stop( 'dyn_dmp')
247      !
248   END SUBROUTINE dyn_dmp
249
250   !!======================================================================
251END MODULE dyndmp
Note: See TracBrowser for help on using the repository browser.