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.
ldfdyn_smag.F90 in branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 4 years ago

The Dr Hook changes from my perl code.

File size: 12.6 KB
Line 
1MODULE ldfdyn_smag
2   !!======================================================================
3   !!                     ***  MODULE  ldftrasmag  ***
4   !! Ocean physics:  variable eddy induced velocity coefficients
5   !!======================================================================
6#if   defined key_dynldf_smag   &&   defined key_dynldf_c3d
7   !!----------------------------------------------------------------------
8   !!   'key_dynldf_smag'      and           smagorinsky  diffusivity
9   !!   'key_dynldf_c3d'                    3D tracer lateral  mixing coef.
10   !!----------------------------------------------------------------------
11   !!   ldf_eiv      : compute the eddy induced velocity coefficients
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE sbc_oce         ! surface boundary condition: ocean
17   USE sbcrnf          ! river runoffs
18   USE ldfdyn_oce      ! ocean tracer   lateral physics
19   USE phycst          ! physical constants
20   USE ldfslp          ! iso-neutral slopes
21   USE in_out_manager  ! I/O manager
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE prtctl          ! Print control
24   USE iom
25   USE wrk_nemo
26   USE yomhook, ONLY: lhook, dr_hook
27   USE parkind1, ONLY: jprb, jpim
28
29   IMPLICIT NONE
30   PRIVATE
31
32   !! * Routine accessibility
33   PUBLIC ldf_dyn_smag               ! routine called by step.F90
34   !!----------------------------------------------------------------------
35   !!  OPA 9.0 , LOCEAN-IPSL (2005)
36   !! $Id$
37   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
38   !!----------------------------------------------------------------------
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43
44CONTAINS
45
46
47
48
49
50   !!----------------------------------------------------------------------
51   !!                        ***  ldfdyn_smag.F90  ***
52   !!----------------------------------------------------------------------
53
54   !!----------------------------------------------------------------------
55   !!  OPA 9.0 , LOCEAN-IPSL (2005)
56   !! $Id$
57   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
58   !!----------------------------------------------------------------------
59
60   !!----------------------------------------------------------------------
61   !!   'key_dynldf_smag'             3D lateral eddy viscosity coefficients
62   !!----------------------------------------------------------------------
63
64   SUBROUTINE ldf_dyn_smag( kt )
65      !!----------------------------------------------------------------------
66      !!                  ***  ROUTINE ldf_dyn_smag  ***
67      !!                   
68      !! ** Purpose :   initializations of the horizontal ocean physics
69      !!
70      !! ** Method  :   3D eddy viscosity coef.
71      !!    M.Griffies, R.Hallberg AMS, 2000
72      !! for laplacian:
73      !!   Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=3-4
74      !! for bilaplacian:
75      !!   Bsmag=Asmag*dx*dy/8
76      !!   D^2=(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates
77      !!  in general case du/dx ==> e2 d(u/e2)/dx;  du/dy ==> e1 d(u/e1)/dy;
78      !!                  dv/dx ==> e2 d(v/e2)/dx;  dv/dy ==> e1 d(v/e1)/dy
79      !!
80      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
81      !!                              ahm3, ahm4 never used
82      !!       bilaplacian operator : ahm1, ahm2 never used
83      !!                           :  ahm3, ahm4 defined at U- and V-points
84      !!       explanation of the default is missingi
85      !!  last modified : Maria Luneva, September 2011
86      !!----------------------------------------------------------------------
87      !! * Modules used
88      !! ahm0 here is a background viscosity
89
90      !! * Arguments
91
92      !! * local variables
93
94      INTEGER              :: kt                   ! timestep
95
96      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
97      REAL (wp):: zdeltat,zdeltaf,zdeltau,zdeltav  ! temporary scalars
98      REAL (wp), POINTER, DIMENSION (:,:) ::   zux, zuy , zvx ,zvy, zue1, zue2, zve1, zve2 
99      REAL (wp)::  zcmsmag_1, zcmsmag_2 , zcmsh
100      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
101      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
102      REAL(KIND=jprb)               :: zhook_handle
103
104      CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_DYN_SMAG'
105
106      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
107
108
109
110      !!----------------------------------------------------------------------
111
112      CALL wrk_alloc( jpi,jpj,zux,zuy,zvx,zvy )
113      CALL wrk_alloc( jpi,jpj,zue1,zue2,zve1,zve2 )
114
115
116      IF(  kt == nit000 ) THEN
117
118
119      IF(lwp) WRITE(numout,*)
120      IF(lwp) WRITE(numout,*) 'ldf_dyn_smag : 3D lateral eddy viscosity coefficient'
121      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
122     
123      ENDIF
124     
125      zcmsmag_1 = rn_cmsmag_1
126      zcmsmag_2 = rn_cmsmag_2
127      zcmsh     = rn_cmsh
128
129   
130
131
132      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operators
133      ! =================                       whatever its orientation is)
134      IF( ln_dynldf_lap ) THEN
135         ! define ahm1 and ahm2 at the right grid point position
136         ! (USER: modify ahm1 and ahm2 following your desiderata)
137         
138         DO jk=1,jpk
139           zue2(:,:)=un(:,:,jk)/e2u(:,:)
140           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
141           zue1(:,:)=un(:,:,jk)/e1u(:,:)
142           zve2(:,:)=vn(:,:,jk)/e2v(:,:)
143
144 
145           DO jj=2,jpj
146            DO ji=2,jpi
147            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) * zcmsh
148            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk) * zcmsh
149            ENDDO
150           ENDDO
151
152           DO jj=1,jpjm1
153            DO ji=1,jpim1
154            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
155            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
156            ENDDO
157           ENDDO
158             
159          DO jj=2,jpjm1
160           DO ji=2,jpim1
161
162            zdeltat=2._wp /(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2))
163            zdeltaf=2._wp /(e1f(ji,jj)**(-2)+e2f(ji,jj)**(-2))   
164            ahm1(ji,jj,jk)=(zcmsmag_1/rpi)**2*zdeltat*                                         &
165                            sqrt( (zux(ji,jj)-zvy(ji,jj))**2+                                    &
166                     0.0625_wp*(zuy(ji,jj)+zuy(ji,jj-1)+zuy(ji-1,jj)+zuy(ji-1,jj-1)+             &
167                            zvx(ji,jj)+zvx(ji,jj-1)+zvx(ji-1,jj)+zvx(ji-1,jj-1))**2)
168
169            ahm2(ji,jj,jk)=(zcmsmag_1/rpi)**2*zdeltaf*                                         &
170                            sqrt( (zuy(ji,jj)+zvx(ji,jj))**2+                                    &
171                     0.0625_wp*(zux(ji,jj)+zux(ji,jj+1)+zux(ji+1,jj)+zux(ji+1,jj+1)-             &
172                             zvy(ji,jj)-zvy(ji,jj+1)-zvy(ji+1,jj)-zvy(ji+1,jj+1))**2)
173
174            ahm1(ji,jj,jk)=MAX(ahm1(ji,jj,jk),rn_ahm_0_lap)
175            ahm2(ji,jj,jk)=MAX(ahm2(ji,jj,jk),rn_ahm_0_lap)
176
177! stability criteria or upper limit set from namelist
178            ahm1(ji,jj,jk)=MIN(ahm1(ji,jj,jk),zdeltat / (16_wp*rdt),rn_ahm_m_lap)
179            ahm2(ji,jj,jk)=MIN(ahm2(ji,jj,jk),zdeltaf / (16_wp*rdt),rn_ahm_m_lap)
180
181            ENDDO
182           ENDDO
183
184         ENDDO ! jpk
185            ahm1(:,:,jpk) = ahm1(:,:,jpkm1)
186            ahm2(:,:,jpk) = ahm2(:,:,jpkm1)
187            IF(lwp.and.kt==nit000) WRITE(numout,'(36x," ahm ", 7x)')
188            DO jk = 1, jpk
189
190               IF(lwp.and.kt==nit000) WRITE(numout,'(30x,E10.2,8x,i3)') ahm1(jpi/2,jpj/2,jk), jk
191            END DO
192      CALL lbc_lnk( ahm1, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
193      CALL lbc_lnk( ahm2, 'F', 1. )   ! Lateral boundary conditions on ( ahtt )
194
195      ENDIF    ! ln_dynldf_lap
196     
197
198
199      ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
200      ! ================================  whatever its orientation is)
201      ! Here: ahm is proportional to the cube of the maximum of the grid spacing
202      !       in the to horizontal direction
203
204      IF( ln_dynldf_bilap ) THEN
205         DO jk=1,jpk
206           zue2(:,:) = un(:,:,jk)/e2u(:,:)
207           zve1(:,:) = vn(:,:,jk)/e1v(:,:)
208           zue1(:,:) = un(:,:,jk)/e1u(:,:)
209           zve2(:,:) = vn(:,:,jk)/e2v(:,:)
210
211
212           DO jj=2,jpj
213            DO ji=2,jpi
214            zux(ji,jj) = (zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
215            zvy(ji,jj) = (zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
216            ENDDO
217           ENDDO
218
219           DO jj=1,jpjm1
220            DO ji=1,jpim1
221            zuy(ji,jj) = (zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
222            zvx(ji,jj) = (zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
223            ENDDO
224           ENDDO
225
226
227          DO jj=2,jpjm1
228           DO ji=2,jpim1
229            zdeltau = 2._wp/(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2))
230            zdeltav = 2._wp/(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2))
231
232            ahm3(ji,jj,jk) = -(zcmsmag_2/rpi)**2/8.0_wp*zdeltau**2*                          &
233
234                         sqrt(0.25_wp*(zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj))**2+    &
235                              0.25_wp*(zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1))**2)
236
237            ahm4(ji,jj,jk) = -(zcmsmag_2/rpi)**2/8.0_wp*zdeltav**2*                           &
238
239                         sqrt(0.25_wp*(zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1))**2+    &
240                              0.25_wp*(zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj))**2)
241
242            ahm3(ji,jj,jk) = MIN (rn_ahm_0_blp , ahm3(ji,jj,jk) )
243            ahm4(ji,jj,jk) = MIN (rn_ahm_0_blp , ahm4(ji,jj,jk) ) 
244
245! stability criteria or upper limit set in namelist
246
247            ahm3(ji,jj,jk) = MAX( ahm3(ji,jj,jk),-zdeltau**2/( 128._wp*rdt ),rn_ahm_m_blp )
248            ahm4(ji,jj,jk) = MAX( ahm4(ji,jj,jk),-zdeltav**2/( 128._wp*rdt ),rn_ahm_m_blp )
249           
250
251            ENDDO
252           ENDDO
253
254         ENDDO
255            ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
256            ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
257
258      DO jk = 1, jpk
259      IF(  kt == nit000 ) THEN
260
261               IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(jpi/2,jpj/2,jk), jk
262      ENDIF   
263      END DO
264      CALL lbc_lnk( ahm3, 'U', 1. )   ! Lateral boundary conditions
265      CALL lbc_lnk( ahm4, 'V', 1. )
266   ENDIF
267
268      CALL wrk_dealloc( jpi,jpj,zux,zuy,zvx,zvy )
269      CALL wrk_dealloc( jpi,jpj,zue1,zue2,zve1,zve2 )
270   !    zumax = MAXVAL( ABS( ahm3(:,:,:) ) )                ! slower than the following loop on NEC SX5
271      zdeltat = 0._wp
272   If(ln_dynldf_lap)THEN
273      DO jk = 1, jpk
274         DO jj = 1, jpj
275            DO ji = 1, jpi
276               zdeltat = MAX(zdeltat,ABS(ahm1(ji,jj,jk)),ABS(ahm2(ji,jj,jk)) )
277          END DO
278        END DO
279      END DO
280      IF( lk_mpp )   CALL mpp_max( zdeltat )                 ! max over the global domain
281      !
282      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,'dynlap:  abs(ahm) max: ', zdeltat
283    ENDIF
284    If(ln_dynldf_bilap)THEN
285     zdeltat = 0._wp
286      DO jk = 1, jpk
287         DO jj = 1, jpj
288            DO ji = 1, jpi
289               zdeltat = MAX(zdeltat,ABS(ahm3(ji,jj,jk)),ABS(ahm3(ji,jj,jk)) )
290          END DO
291        END DO
292      END DO
293      IF( lk_mpp )   CALL mpp_max( zdeltat )                 ! max over the global domain
294      !
295      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,'dyn_bilap abs(ahm) max: ', zdeltat
296      !
297   ENDIF
298      !
299
300      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
301END SUBROUTINE ldf_dyn_smag
302#else
303   !!----------------------------------------------------------------------
304   !!   Default option                                         Dummy module
305   !!----------------------------------------------------------------------
306CONTAINS
307   SUBROUTINE ldf_dyn_smag( kt )       ! Empty routine
308   USE yomhook, ONLY: lhook, dr_hook
309   USE parkind1, ONLY: jprb, jpim
310
311   IMPLICIT NONE
312      INTEGER :: kt                    ! timestep   
313      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
314      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
315      REAL(KIND=jprb)               :: zhook_handle
316
317      CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_DYN_SMAG'
318
319      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
320
321      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
322      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
323   END SUBROUTINE ldf_dyn_smag
324#endif
325
326   END MODULE ldfdyn_smag
327
Note: See TracBrowser for help on using the repository browser.