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 trunk/NEMOGCM/NEMO/OPA_SRC/C1D – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

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