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_c2d.h90 in trunk/NEMO/OPA_SRC/LDF – NEMO

source: trunk/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 10.9 KB
Line 
1   !!----------------------------------------------------------------------
2   !!                      ***  ldfdyn_c2d.h90  ***
3   !!----------------------------------------------------------------------
4   !!   ldf_dyn_c2d  : set the lateral viscosity coefficients
5   !!   ldf_dyn_c2d_orca : specific case for orca r2 and r4
6   !!----------------------------------------------------------------------
7
8   !!----------------------------------------------------------------------
9   !!   OPA 9.0 , LOCEAN-IPSL (2005)
10   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90,v 1.6 2007/06/29 17:01:51 opalod Exp $
11   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
12   !!----------------------------------------------------------------------
13
14   SUBROUTINE ldf_dyn_c2d( ld_print )
15      !!----------------------------------------------------------------------
16      !!                 ***  ROUTINE ldf_dyn_c2d  ***
17      !!                 
18      !! ** Purpose :   initializations of the horizontal ocean physics
19      !!
20      !! ** Method :
21      !!      2D eddy viscosity coefficients ( longitude, latitude )
22      !!
23      !!       harmonic operator   : ahm1 is defined at t-point
24      !!                             ahm2 is defined at f-point
25      !!           + isopycnal     : ahm3 is defined at u-point
26      !!           or geopotential   ahm4 is defined at v-point
27      !!           iso-model level : ahm3, ahm4 not used
28      !!
29      !!       biharmonic operator : ahm1 is defined at u-point
30      !!                             ahm2 is defined at v-point
31      !!                           : ahm3, ahm4 not used
32      !!
33      !!----------------------------------------------------------------------
34      !! * Arguments
35      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
36
37      !! * Local variables
38      REAL(wp) ::   za00, zdx_max
39      !!----------------------------------------------------------------------
40
41      IF(lwp) WRITE(numout,*)
42      IF(lwp) WRITE(numout,*) 'ldf_dyn_c2d : 2d lateral eddy viscosity coefficient'
43      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
44      IF(lwp) WRITE(numout,*)
45
46      ! harmonic operator (ahm1, ahm2) : ( T- and F- points) (used for laplacian operators
47      ! ===============================                       whatever its orientation is)
48      IF( ln_dynldf_lap ) THEN
49         ! define ahm1 and ahm2 at the right grid point position
50         ! (USER: modify ahm1 and ahm2 following your desiderata)
51
52         zdx_max = MAXVAL( e1t(:,:) )
53         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
54
55         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
56         IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
57         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
58
59         za00 = ahm0 / zdx_max
60         ahm1(:,:) = za00 * e1t(:,:)
61         ahm2(:,:) = za00 * e1f(:,:)
62
63         IF( ln_dynldf_iso ) THEN
64            IF(lwp) WRITE(numout,*) '              Caution, as implemented now, the isopycnal part of momentum'
65            IF(lwp) WRITE(numout,*) '                 mixing use aht0 as eddy viscosity coefficient. Thus, it is'
66            IF(lwp) WRITE(numout,*) '                 uniform and you must be sure that your ahm is greater than'
67            IF(lwp) WRITE(numout,*) '                 aht0 everywhere in the model domain.'
68         ENDIF
69
70         ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahm1 ahm2)
71         ! ==============================================
72         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   CALL ldf_dyn_c2d_orca( ld_print )
73
74         ! Control print
75         IF( lwp .AND. ld_print ) THEN
76            WRITE(numout,*)
77            WRITE(numout,*) 'inildf: 2D ahm1 array'
78            CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
79            WRITE(numout,*)
80            WRITE(numout,*) 'inildf: 2D ahm2 array'
81            CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
82         ENDIF
83      ENDIF
84
85      ! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator
86      ! =================================                      whatever its orientation is)
87      IF( ln_dynldf_bilap ) THEN
88         ! (USER: modify ahm3 and ahm4 following your desiderata)
89         ! Here: ahm is proportional to the cube of the maximum of the gridspacing
90         !       in the to horizontal direction
91
92         zdx_max = MAXVAL( e1u(:,:) )
93         IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain
94
95         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
96         IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...'
97         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0
98
99         za00 = ahm0 / ( zdx_max * zdx_max * zdx_max )
100         ahm3(:,:) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:)
101         ahm4(:,:) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:)
102
103         ! Control print
104         IF( lwp .AND. ld_print ) THEN
105            WRITE(numout,*)
106            WRITE(numout,*) 'inildf: ahm3 array'
107            CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
108            WRITE(numout,*)
109            WRITE(numout,*) 'inildf: ahm4 array'
110            CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
111         ENDIF
112      ENDIF
113
114
115   END SUBROUTINE ldf_dyn_c2d
116
117
118   SUBROUTINE ldf_dyn_c2d_orca( ld_print )
119      !!----------------------------------------------------------------------
120      !!                 ***  ROUTINE ldf_dyn_c2d  ***
121      !!
122      !!                   **** W A R N I N G ****
123      !!
124      !!                ORCA R2 and R4 configurations
125      !!                 
126      !!                   **** W A R N I N G ****
127      !!                 
128      !! ** Purpose :   initializations of the lateral viscosity for orca R2
129      !!
130      !! ** Method  :   blah blah blah...
131      !!
132      !!----------------------------------------------------------------------
133      !! * Modules used
134      USE ldftra_oce, ONLY : aht0
135
136      !! * Arguments
137      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
138
139      !! * Local variables
140      INTEGER ::   ji, jj, jn      ! dummy loop indices
141      INTEGER ::   inum            ! temporary logical unit
142      INTEGER ::   iim, ijm
143      INTEGER ::   ifreq, il1, il2, ij, ii
144      INTEGER, DIMENSION(jpidta,jpidta) ::   idata
145      INTEGER, DIMENSION(jpi   ,jpj   ) ::   icof
146
147      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
148
149      CHARACTER (len=15) ::   clexp
150      !!----------------------------------------------------------------------
151
152      IF(lwp) WRITE(numout,*)
153      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
154      IF(lwp) WRITE(numout,*) '~~~~~~  --'
155      IF(lwp) WRITE(numout,*)
156      IF(lwp) WRITE(numout,*) '        orca ocean model'
157      IF(lwp) WRITE(numout,*)
158
159#if defined key_antarctic
160#     include "ldfdyn_antarctic.h90"
161#elif defined key_arctic
162#     include "ldfdyn_arctic.h90"
163#else
164      ! Read 2d integer array to specify western boundary increase in the
165      ! ===================== equatorial strip (20N-20S) defined at t-points
166
167      CALL ctlopn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
168         &           1, numout, lwp, 1 )
169      REWIND inum
170      READ(inum,9101) clexp, iim, ijm
171      READ(inum,'(/)')
172      ifreq = 40
173      il1 = 1
174      DO jn = 1, jpidta/ifreq+1
175         READ(inum,'(/)')
176         il2 = MIN( jpidta, il1+ifreq-1 )
177         READ(inum,9201) ( ii, ji = il1, il2, 5 )
178         READ(inum,'(/)')
179         DO jj = jpjdta, 1, -1
180            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
181         END DO
182         il1 = il1 + ifreq
183      END DO
184     
185      DO jj = 1, nlcj
186         DO ji = 1, nlci
187            icof(ji,jj) = idata( mig(ji), mjg(jj) )
188         END DO
189      END DO
190      DO jj = nlcj+1, jpj
191         DO ji = 1, nlci
192            icof(ji,jj) = icof(ji,nlcj)
193         END DO
194      END DO
195      DO jj = 1, jpj
196         DO ji = nlci+1, jpi
197            icof(ji,jj) = icof(nlci,jj)
198         END DO
199      END DO
200
201 9101 FORMAT(1x,a15,2i8)
202 9201 FORMAT(3x,13(i3,12x))
203 9202 FORMAT(i3,41i3)
204
205
206      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
207      ! =================
208      ! define ahm1 and ahm2 at the right grid point position
209      ! (USER: modify ahm1 and ahm2 following your desiderata)
210     
211     
212      ! Decrease ahm to zahmeq m2/s in the tropics
213      ! (from 90 to 20 degre: ahm = constant
214      ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
215      ! from  2.5 to  0 degre: ahm = constant
216      ! symmetric in the south hemisphere)
217
218      zahmeq = aht0
219     
220      DO jj = 1, jpj
221         DO ji = 1, jpi
222            IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
223               ahm2(ji,jj) =  ahm0
224            ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
225               ahm2(ji,jj) =  zahmeq
226            ELSE
227               ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
228                  * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
229            ENDIF
230            IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
231               ahm1(ji,jj) =  ahm0
232            ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
233               ahm1(ji,jj) =  zahmeq
234            ELSE
235               ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
236                  * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
237            ENDIF
238         END DO
239      END DO
240
241      ! increase along western boundaries of equatorial strip
242      ! t-point
243      DO jj = 1, jpjm1
244         DO ji = 1, jpim1
245            zcoft = FLOAT( icof(ji,jj) ) / 100.
246            ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
247         END DO
248      END DO
249      ! f-point
250      icof(:,:) = icof(:,:) * tmask(:,:,1)
251      DO jj = 1, jpjm1
252         DO ji = 1, jpim1
253            zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
254            IF( zmsk == 0. ) THEN
255               zcoff = 1.
256            ELSE
257               zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
258                     / (zmsk * 100.)
259            ENDIF
260            ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
261         END DO
262      END DO
263#endif
264     
265      ! Lateral boundary conditions on ( ahm1, ahm2 )
266      !                                ==============
267      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
268      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
269
270      ! Control print
271      IF( lwp .AND. ld_print ) THEN
272         WRITE(numout,*)
273         WRITE(numout,*) 'inildf: 2D ahm1 array'
274         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
275         WRITE(numout,*)
276         WRITE(numout,*) 'inildf: 2D ahm2 array'
277         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
278      ENDIF
279
280   END SUBROUTINE ldf_dyn_c2d_orca
Note: See TracBrowser for help on using the repository browser.