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

source: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90 @ 7424

Last change on this file since 7424 was 2932, checked in by hliu, 13 years ago
File size: 10.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
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: ldf_tra_smag.F90 1482 2010-06-13 15:28:06Z  $
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: ldfdyn_c3d.h90 1581 2009-08-05 14:53:12Z smasson $
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      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
90      USE wrk_nemo, ONLY:   zux => wrk_2d_1  , zuy => wrk_2d_2  , zvx => wrk_2d_3  , zvy => wrk_2d_4   ! 2D workspace
91      USE wrk_nemo, ONLY:   zue1 => wrk_2d_5 , zue2 => wrk_2d_6 , zve1 => wrk_2d_7 , zve2 => wrk_2d_8   ! 2D workspac
92      !! * local variables
93
94      INTEGER              :: kt                   ! timestep
95
96      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
97      REAL (wp):: rdeltat,rdeltaf,rdeltau,rdeltav  ! temporary scalars
98     
99
100      !!----------------------------------------------------------------------
101
102      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
103         CALL ctl_stop('dyn_ldf_smag : requested workspace arrays unavailable')   ;   RETURN
104      ENDIF
105
106
107      IF(  kt == nit000 ) THEN
108
109
110      IF(lwp) WRITE(numout,*)
111      IF(lwp) WRITE(numout,*) 'ldf_dyn_smag : 3D lateral eddy viscosity coefficient'
112      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
113     
114      ENDIF
115   
116
117
118      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operators
119      ! =================                       whatever its orientation is)
120      IF( ln_dynldf_lap ) THEN
121         ! define ahm1 and ahm2 at the right grid point position
122         ! (USER: modify ahm1 and ahm2 following your desiderata)
123         
124         DO jk=1,jpk
125           zue2(:,:)=un(:,:,jk)/e2u(:,:)
126           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
127           zue1(:,:)=un(:,:,jk)/e1u(:,:)
128           zve2(:,:)=vn(:,:,jk)/e2v(:,:)
129
130 
131           DO jj=2,jpj
132            DO ji=2,jpi
133            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
134            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
135            ENDDO
136           ENDDO
137
138           DO jj=1,jpjm1
139            DO ji=1,jpi
140            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
141            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
142            ENDDO
143           ENDDO
144             
145          DO jj=2,jpjm1
146           DO ji=2,jpim1
147            rdeltat=2./(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2))
148            rdeltaf=2./(e1f(ji,jj)**(-2)+e2f(ji,jj)**(-2))   
149            ahm1(ji,jj,jk)=(rn_cmsmag_1/3.14)**2*rdeltat*                                     &
150                            sqrt( (zux(ji,jj)-zvy(ji,jj))**2+                                 &
151                     0.0625*(zuy(ji,jj)+zuy(ji,jj-1)+zuy(ji-1,jj)+zuy(ji-1,jj-1)+             &
152                            zvx(ji,jj)+zvx(ji,jj-1)+zvx(ji-1,jj)+zvx(ji-1,jj-1))**2)
153
154            ahm2(ji,jj,jk)=(rn_cmsmag_1/3.14)**2*rdeltaf*                                     &
155                            sqrt( (zuy(ji,jj)+zvx(ji,jj))**2+                                 &
156                     0.0625*(zux(ji,jj)+zux(ji,jj+1)+zux(ji+1,jj)+zux(ji+1,jj+1)-             &
157                             zvy(ji,jj)-zvy(ji,jj+1)-zvy(ji+1,jj)-zvy(ji+1,jj+1))**2)
158
159            ahm1(ji,jj,jk)=MAX(ahm1(ji,jj,jk),rn_ahm_0_lap)
160            ahm2(ji,jj,jk)=MAX(ahm2(ji,jj,jk),rn_ahm_0_lap)
161
162! stability criteria
163            ahm1(ji,jj,jk)=MIN(ahm1(ji,jj,jk),rdeltat**2/(16*rdt))
164            ahm2(ji,jj,jk)=MIN(ahm2(ji,jj,jk),rdeltaf**2/(16*rdt))
165
166            ENDDO
167           ENDDO
168
169         ENDDO ! jpk
170            ahm1(:,:,jpk) = ahm1(:,:,jpkm1)
171            ahm2(:,:,jpk) = ahm2(:,:,jpkm1)
172            IF(lwp.and.kt==nit000) WRITE(numout,'(36x," ahm ", 7x)')
173            DO jk = 1, jpk
174
175               IF(lwp.and.kt==nit000) WRITE(numout,'(30x,E10.2,8x,i3)') ahm1(jpi/2,jpj/2,jk), jk
176            END DO
177      CALL lbc_lnk( ahm1, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
178      CALL lbc_lnk( ahm2, 'F', 1. )   ! Lateral boundary conditions on ( ahtt )
179
180      ENDIF    ! ln_dynldf
181     
182
183
184      ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
185      ! ================================  whatever its orientation is)
186      ! (USER: modify ahm3 and ahm4 following your desiderata)
187      ! Here: ahm is proportional to the cube of the maximum of the gridspacing
188      !       in the to horizontal direction
189
190      IF( ln_dynldf_bilap ) THEN
191         DO jk=1,jpk
192           zue2(:,:)=un(:,:,jk)/e2u(:,:)
193           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
194           zue1(:,:)=un(:,:,jk)/e1u(:,:)
195           zve2(:,:)=vn(:,:,jk)/e2v(:,:)
196
197
198           DO jj=2,jpj
199            DO ji=2,jpi
200            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
201            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
202            ENDDO
203           ENDDO
204
205           DO jj=1,jpjm1
206            DO ji=1,jpim1
207            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
208            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
209            ENDDO
210           ENDDO
211
212
213          DO jj=2,jpjm1
214           DO ji=2,jpim1
215            rdeltau=2./(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2))
216            rdeltav=2./(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2))
217
218             ahm3(ji,jj,jk)=MIN(rn_ahm_0_blp,(rn_cmsmag_2/3.14)**2/8*rdeltau**2*              &
219
220                         sqrt(0.25*(zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj))**2+    &
221                              0.25*(zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1))**2))
222
223            ahm4(ji,jj,jk)=MIN(rn_ahm_0_blp ,(rn_cmsmag_2/3.14)**2/8*rdeltav**2*              &
224
225                         sqrt(0.25*(zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1))**2+    &
226                              0.25*(zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj))**2))
227! stability criteria
228            ahm3(ji,jj,jk)=MAX(ahm3(ji,jj,jk),-rdeltau**2/(16*rdt))
229            ahm4(ji,jj,jk)=MAX(ahm4(ji,jj,jk),-rdeltav**2/(16*rdt))
230           
231
232            ENDDO
233           ENDDO
234
235         ENDDO
236            ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
237            ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
238
239      DO jk = 1, jpk
240      IF(  kt == nit000 ) THEN
241
242               IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(jpi/2,jpj/2,jk), jk
243      ENDIF   
244      END DO
245      CALL lbc_lnk( ahm3, 'U', 1. )   ! Lateral boundary conditions
246      CALL lbc_lnk( ahm4, 'V', 1. )
247   ENDIF
248
249   IF( wrk_not_released(2, 1,2,3,4,5,6,7,8)) CALL ctl_stop('dyn_ldf_smag: failed to release workspace arrays')
250END SUBROUTINE ldf_dyn_smag
251#else
252   !!----------------------------------------------------------------------
253   !!   Default option                                         Dummy module
254   !!----------------------------------------------------------------------
255CONTAINS
256   SUBROUTINE ldf_dyn_smag( kt )       ! Empty routine
257      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
258   END SUBROUTINE ldf_dyn_smag
259#endif
260
261   END MODULE ldfdyn_smag
262
Note: See TracBrowser for help on using the repository browser.