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 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

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