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_package/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90 @ 13190

Last change on this file since 13190 was 11101, checked in by frrh, 5 years ago

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