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.
ldftra_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/ldftra_smag.F90 @ 11101

Last change on this file since 11101 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: 9.1 KB
Line 
1MODULE ldftra_smag
2   !!======================================================================
3   !!                     ***  MODULE  ldftrasmag  ***
4   !! Ocean physics:  variable eddy induced velocity coefficients
5   !!======================================================================
6#if   defined key_traldf_smag   &&   defined key_traldf_c3d
7   !!----------------------------------------------------------------------
8   !!   'key_traldf_smag'      and           smagorinsky  diffusivity
9   !!   'key_traldf_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 ldftra_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_tra_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  !!                        ***  ldf_tra_smag.F90  ***
49  !!----------------------------------------------------------------------
50
51
52   SUBROUTINE ldf_tra_smag( kt )
53      !!----------------------------------------------------------------------
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE ldf_tra_smag  ***
56      !!                   
57      !! ** Purpose :   initializations of the horizontal ocean physics
58      !!
59      !! ** Method  :   3D eddy viscosity coef.
60      !!    M.Griffies, R.Hallberg AMS, 2000
61      !!  for laplacian:
62      !!   Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=1 for tracers, C=3-4 for viscosity
63      !!   D^2= rm_smsh*(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates
64      !!  IF rm_smsh = 0 , only shear is used, recommended for tidal flows
65      !!  in general case du/dx ==> e2 d(u/e2)/dx;  du/dy ==> e1 d(u/e1)/dy;
66      !!                  dv/dx ==> e2 d(v/e2)/dx;  dv/dy ==> e1 d(v/e1)/dy
67      !!  for bilaplacian: now this option is deleted as unstable or non-conservative
68      !!   - delta{Bsmag (delta(T)} = -Bsmag* delta{delta(T)} - delta(Bsmag)*delta( T )
69      !!  second term is of arbitrary sign on the edge of fronts and can induce instability
70      !!   Bsmag=Asmag*dx*dy/8
71      !! 
72      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
73      !!                              ahm3, ahm4 never used
74      !!       bilaplacian operator : ahm1, ahm2 never used
75      !!                           :  ahm3, ahm4 defined at U- and V-points
76      !!       ??? explanation of the default is missing
77      !!  last modified : Maria Luneva, October 2012
78      !!----------------------------------------------------------------------
79      !!
80      !!----------------------------------------------------------------------
81      !! * Modules used
82      USE ioipsl
83      REAL ( wp), POINTER , DIMENSION (:,:) :: zux, zvx , zuy , zvy 
84      REAL ( wp), POINTER , DIMENSION (:,:) :: zue1, zue2 , zve1 , zve2 
85      INTEGER, INTENT( in )                 ::   kt                             ! ocean time-step inedx
86      !! * Arguments
87      INTEGER                               :: ji,jj,jk
88
89      REAL (wp)                             ::     zdeltau, zdeltav, zhsmag ,zsmsh    ! temporary scalars
90     
91      CALL wrk_alloc (jpi,jpj,zux, zvx , zuy , zvy )
92      CALL wrk_alloc (jpi,jpj,zue1, zue2 , zve1 , zve2 )
93      !!----------------------------------------------------------------------
94      IF(  kt == nit000 ) THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity '
97         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  '
98         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*)
101         IF(lwp .AND. lflush) CALL flush(numout)
102      ENDIF
103
104      zhsmag = rn_chsmag
105      zsmsh  = rn_smsh 
106      zux(:,:)=0._wp ; zuy(:,:)=0._wp ; zvx(:,:)=0._wp ; zvy(:,:)=0._wp
107
108      ! -------------------
109      ahtt(:,:,:) = rn_aht_0
110       IF( ln_traldf_bilap ) THEN
111        IF( lwp .AND. kt == nit000) THEN
112          WRITE(numout,* )'ldf_tra_smag :no bilaplacian Smagorinsky diffusivity'
113          WRITE(numout,* )'ldf_tra_smag :bilaplacian diffusivity set to constant' 
114          IF(lflush) CALL flush(numout)
115        ENDIF
116       ENDIF
117
118
119
120      ! harmonic operator   (U-, V-, W-points)
121      ! -----------------
122
123      ahtu(:,:,:) = rn_aht_0                  ! set ahtu , ahtv at u- and v-points,
124      ahtv(:,:,:) = rn_aht_0                  ! and ahtw at w-point
125      ahtw(:,:,:) = rn_aht_0                  ! (here example: no space variation)
126     
127      IF( ln_traldf_lap ) THEN
128
129         DO jk=1,jpk
130
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) * zsmsh 
140            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk) * zsmsh 
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
152          DO jj=2,jpjm1
153           DO ji=2,jpim1
154           zdeltau=2._wp/( e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2) )
155           zdeltav=2._wp/( e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2) )
156
157           ahtu(ji,jj,jk)=MAX( rn_aht_0 , (zhsmag/rpi)**2*zdeltau*                                &
158                          SQRT(0.25_wp*( zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj) )**2+    &
159                               0.25_wp*( zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1)  )**2) )
160
161           ahtv(ji,jj,jk)=MAX( rn_aht_0 ,  (zhsmag/rpi)**2*zdeltav*                               &
162                          SQRT(0.25_wp*( zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1) )**2+    &
163                               0.25_wp*( zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj)  )**2) )
164
165
166         !!! stability criteria: aht<delta**2/(4*dt)   dt=2*rdt , positiveness require aht<delta**2/(8*dt)
167             ahtu(ji,jj,jk)=MIN(ahtu(ji,jj,jk),zdeltau/(16*rdt) ,rn_aht_m)
168             ahtv(ji,jj,jk)=MIN(ahtv(ji,jj,jk),zdeltav/(16*rdt) ,rn_aht_m)
169         ! so...
170
171
172            ENDDO
173           ENDDO
174         ENDDO
175        ENDIF
176            ahtu(:,:,jpk) = ahtu(:,:,jpkm1)
177            ahtv(:,:,jpk) = ahtv(:,:,jpkm1)
178        CALL lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions
179        CALL lbc_lnk( ahtv, 'V', 1. )
180
181IF(  kt == nit000 ) THEN
182
183      IF(lwp ) THEN                    ! Control print
184         WRITE(numout,*)
185         WRITE(numout,*) 'inildf: ahtu at k = 1'
186         CALL prihre( ahtu(:,:,1), jpi, jpj, 1, jpi, 1,   &
187            &                                1, jpj, 1, 1.e-1, numout )
188         WRITE(numout,*)
189         WRITE(numout,*) 'inildf: ahtv at k = 1'
190         CALL prihre( ahtv(:,:,1), jpi, jpj, 1, jpi, 1,   &
191            &                                1, jpj, 1, 1.e-1, numout )
192         WRITE(numout,*)
193         WRITE(numout,*) 'inildf: ahtw at k = 1'
194         CALL prihre( ahtw(:,:,1), jpi, jpj, 1, jpi, 1,   &
195            &                                1, jpj, 1, 1.e-1, numout )
196         IF(lflush) CALL flush(numout)
197      ENDIF
198ENDIF
199
200      CALL wrk_dealloc ( jpi,jpj,zux, zvx , zuy , zvy     )
201      CALL wrk_dealloc ( jpi,jpj,zue1, zue2 , zve1 , zve2 )
202
203
204END SUBROUTINE ldf_tra_smag
205#else
206   !!----------------------------------------------------------------------
207   !!   Default option                                         Dummy module
208   !!----------------------------------------------------------------------
209CONTAINS
210   SUBROUTINE ldf_tra_smag( kt )       ! Empty routine
211   IMPLICIT NONE
212        INTEGER, INTENT( in ) ::   kt  ! ocean time-step inedx
213        WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
214   END SUBROUTINE ldf_tra_smag
215#endif
216
217END MODULE ldftra_smag
Note: See TracBrowser for help on using the repository browser.