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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 2698

Last change on this file since 2698 was 2698, checked in by rblod, 12 years ago

Suppress unmatched END IF

  • Property svn:keywords set to Id
File size: 17.1 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   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
10   !! $Id$
11   !! Software governed by the CeCILL licence (NEMOGCM/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 : ahm3 is defined at u-point
30      !!                             ahm4 is defined at v-point
31      !!                           : ahm1, ahm2 not used
32      !!
33      !!----------------------------------------------------------------------
34      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
35      !
36      INTEGER  ::   ji, jj
37      REAL(wp) ::   za00, zd_max, zetmax, zeumax, zefmax, zevmax
38      !!----------------------------------------------------------------------
39
40      IF(lwp) WRITE(numout,*)
41      IF(lwp) WRITE(numout,*) 'ldf_dyn_c2d : 2d lateral eddy viscosity coefficient'
42      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
43
44      ! harmonic operator (ahm1, ahm2) : ( T- and F- points) (used for laplacian operators
45      ! ===============================                       whatever its orientation is)
46      IF( ln_dynldf_lap ) THEN
47         ! define ahm1 and ahm2 at the right grid point position
48         ! (USER: modify ahm1 and ahm2 following your desiderata)
49
50         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
51         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
52
53         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
54         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
55
56         za00 = ahm0 / zd_max
57         DO jj = 1, jpj
58            DO ji = 1, jpi
59               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
60               zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
61               ahm1(ji,jj) = za00 * zetmax
62               ahm2(ji,jj) = za00 * zefmax
63            END DO
64         END DO
65
66         IF( ln_dynldf_iso ) THEN
67            IF(lwp) WRITE(numout,*) '              Caution, as implemented now, the isopycnal part of momentum'
68            IF(lwp) WRITE(numout,*) '                 mixing use aht0 as eddy viscosity coefficient. Thus, it is'
69            IF(lwp) WRITE(numout,*) '                 uniform and you must be sure that your ahm is greater than'
70            IF(lwp) WRITE(numout,*) '                 aht0 everywhere in the model domain.'
71         ENDIF
72
73         ! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2)
74         ! ==============================================
75         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   CALL ldf_dyn_c2d_orca( ld_print )
76         IF( cp_cfg == "orca" .AND.   jp_cfg == 1)                       CALL ldf_dyn_c2d_orca_R1( ld_print )
77
78         ! Control print
79         IF( lwp .AND. ld_print ) THEN
80            WRITE(numout,*)
81            WRITE(numout,*) 'inildf: 2D ahm1 array'
82            CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
83            WRITE(numout,*)
84            WRITE(numout,*) 'inildf: 2D ahm2 array'
85            CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
86         ENDIF
87      ENDIF
88
89      ! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator
90      ! =================================                      whatever its orientation is)
91      IF( ln_dynldf_bilap ) THEN
92         ! (USER: modify ahm3 and ahm4 following your desiderata)
93         ! Here: ahm is proportional to the cube of the maximum of the gridspacing
94         !       in the to horizontal direction
95
96         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
97         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
98
99         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
100         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
101
102         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
103         DO jj = 1, jpj
104            DO ji = 1, jpi
105               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
106               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
107               ahm3(ji,jj) = za00 * zeumax * zeumax * zeumax
108               ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax
109            END DO
110         END DO
111
112         ! Control print
113         IF( lwp .AND. ld_print ) THEN
114            WRITE(numout,*)
115            WRITE(numout,*) 'inildf: ahm3 array'
116            CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
117            WRITE(numout,*)
118            WRITE(numout,*) 'inildf: ahm4 array'
119            CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
120         ENDIF
121      ENDIF
122      !
123   END SUBROUTINE ldf_dyn_c2d
124
125
126   SUBROUTINE ldf_dyn_c2d_orca( ld_print )
127      !!----------------------------------------------------------------------
128      !!                 ***  ROUTINE ldf_dyn_c2d  ***
129      !!
130      !!                   **** W A R N I N G ****
131      !!
132      !!                ORCA R2 and R4 configurations
133      !!                 
134      !!                   **** W A R N I N G ****
135      !!                 
136      !! ** Purpose :   initializations of the lateral viscosity for orca R2
137      !!
138      !! ** Method  :   blah blah blah...
139      !!
140      !!----------------------------------------------------------------------
141      USE ldftra_oce, ONLY:   aht0
142      USE wrk_nemo  , ONLY:   iwrk_in_use, iwrk_not_released
143      USE wrk_nemo  , ONLY:   icof => iwrk_2d_1
144      !
145      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
146      !
147      INTEGER  ::   ji, jj, jn   ! dummy loop indices
148      INTEGER  ::   inum, iim, ijm            ! local integers
149      INTEGER  ::   ifreq, il1, il2, ij, ii
150      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
151      CHARACTER (len=15) ::   clexp
152      INTEGER, DIMENSION(jpidta,jpidta) ::   idata
153      !!----------------------------------------------------------------------
154
155      IF( iwrk_in_use(2, 1) )THEN
156         CALL ctl_stop('ldf_dyn_c2d_orca: requested workspace array is unavailable')   ;   RETURN
157      ENDIF
158
159      IF(lwp) WRITE(numout,*)
160      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
161      IF(lwp) WRITE(numout,*) '~~~~~~  --'
162      IF(lwp) WRITE(numout,*) '        orca ocean configuration'
163
164#if defined key_antarctic
165#     include "ldfdyn_antarctic.h90"
166#elif defined key_arctic
167#     include "ldfdyn_arctic.h90"
168#else
169      ! Read 2d integer array to specify western boundary increase in the
170      ! ===================== equatorial strip (20N-20S) defined at t-points
171
172      CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
173      READ(inum,9101) clexp, iim, ijm
174      READ(inum,'(/)')
175      ifreq = 40
176      il1 = 1
177      DO jn = 1, jpidta/ifreq+1
178         READ(inum,'(/)')
179         il2 = MIN( jpidta, il1+ifreq-1 )
180         READ(inum,9201) ( ii, ji = il1, il2, 5 )
181         READ(inum,'(/)')
182         DO jj = jpjdta, 1, -1
183            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
184         END DO
185         il1 = il1 + ifreq
186      END DO
187     
188      DO jj = 1, nlcj
189         DO ji = 1, nlci
190            icof(ji,jj) = idata( mig(ji), mjg(jj) )
191         END DO
192      END DO
193      DO jj = nlcj+1, jpj
194         DO ji = 1, nlci
195            icof(ji,jj) = icof(ji,nlcj)
196         END DO
197      END DO
198      DO jj = 1, jpj
199         DO ji = nlci+1, jpi
200            icof(ji,jj) = icof(nlci,jj)
201         END DO
202      END DO
203
204 9101 FORMAT(1x,a15,2i8)
205 9201 FORMAT(3x,13(i3,12x))
206 9202 FORMAT(i3,41i3)
207
208
209      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
210      ! =================
211      ! define ahm1 and ahm2 at the right grid point position
212      ! (USER: modify ahm1 and ahm2 following your desiderata)
213     
214     
215      ! Decrease ahm to zahmeq m2/s in the tropics
216      ! (from 90 to 20 degre: ahm = constant
217      ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
218      ! from  2.5 to  0 degre: ahm = constant
219      ! symmetric in the south hemisphere)
220
221      zahmeq = aht0
222     
223      DO jj = 1, jpj
224         DO ji = 1, jpi
225            IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
226               ahm2(ji,jj) =  ahm0
227            ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
228               ahm2(ji,jj) =  zahmeq
229            ELSE
230               ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
231                  * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
232            ENDIF
233            IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
234               ahm1(ji,jj) =  ahm0
235            ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
236               ahm1(ji,jj) =  zahmeq
237            ELSE
238               ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
239                  * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
240            ENDIF
241         END DO
242      END DO
243
244      ! increase along western boundaries of equatorial strip
245      ! t-point
246      DO jj = 1, jpjm1
247         DO ji = 1, jpim1
248            zcoft = FLOAT( icof(ji,jj) ) / 100.
249            ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
250         END DO
251      END DO
252      ! f-point
253      icof(:,:) = icof(:,:) * tmask(:,:,1)
254      DO jj = 1, jpjm1
255         DO ji = 1, jpim1   ! NO vector opt.
256            zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
257            IF( zmsk == 0. ) THEN
258               zcoff = 1.
259            ELSE
260               zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
261                     / (zmsk * 100.)
262            ENDIF
263            ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
264         END DO
265      END DO
266#endif
267     
268      ! Lateral boundary conditions on ( ahm1, ahm2 )
269      !                                ==============
270      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
271      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
272
273      ! Control print
274      IF( lwp .AND. ld_print ) THEN
275         WRITE(numout,*)
276         WRITE(numout,*) 'inildf: 2D ahm1 array'
277         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
278         WRITE(numout,*)
279         WRITE(numout,*) 'inildf: 2D ahm2 array'
280         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
281      ENDIF
282
283      IF( iwrk_not_released(2, 1) )   CALL ctl_stop('ldf_dyn_c2d_orca: failed to release workspace array')
284      !
285   END SUBROUTINE ldf_dyn_c2d_orca
286
287
288   SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print )
289      !!----------------------------------------------------------------------
290      !!                 ***  ROUTINE ldf_dyn_c2d  ***
291      !!
292      !!                   **** W A R N I N G ****
293      !!
294      !!                ORCA R1 configuration
295      !!                 
296      !!                   **** W A R N I N G ****
297      !!                 
298      !! ** Purpose :   initializations of the lateral viscosity for orca R1
299      !!
300      !! ** Method  :   blah blah blah...
301      !!
302      !!----------------------------------------------------------------------
303      USE ldftra_oce, ONLY:   aht0
304      USE wrk_nemo  , ONLY:   iwrk_in_use, iwrk_not_released
305      USE wrk_nemo  , ONLY:   icof => iwrk_2d_1
306      !
307      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
308      !
309      INTEGER ::   ji, jj, jn      ! dummy loop indices
310      INTEGER ::   inum            ! temporary logical unit
311      INTEGER ::   iim, ijm
312      INTEGER ::   ifreq, il1, il2, ij, ii
313      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s
314      CHARACTER (len=15) ::   clexp
315      INTEGER, DIMENSION(jpidta,jpidta) ::   idata
316      !!----------------------------------------------------------------------
317
318      IF( iwrk_in_use(2, 1) ) THEN
319         CALL ctl_stop('ldf_dyn_c2d_orca_R1: requested workspace array is unavailable')   ;   RETURN
320      ENDIF
321
322      IF(lwp) WRITE(numout,*)
323      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
324      IF(lwp) WRITE(numout,*) '~~~~~~  --'
325      IF(lwp) WRITE(numout,*) '        orca_r1 configuration'
326
327#if defined key_antarctic
328#     include "ldfdyn_antarctic.h90"
329#elif defined key_arctic
330#     include "ldfdyn_arctic.h90"
331#else
332      ! Read 2d integer array to specify western boundary increase in the
333      ! ===================== equatorial strip (20N-20S) defined at t-points
334
335      CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
336         &           1, numout, lwp )
337      REWIND inum
338      READ(inum,9101) clexp, iim, ijm
339      READ(inum,'(/)')
340      ifreq = 40
341      il1 = 1
342      DO jn = 1, jpidta/ifreq+1
343         READ(inum,'(/)')
344         il2 = MIN( jpidta, il1+ifreq-1 )
345         READ(inum,9201) ( ii, ji = il1, il2, 5 )
346         READ(inum,'(/)')
347         DO jj = jpjdta, 1, -1
348            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
349         END DO
350         il1 = il1 + ifreq
351      END DO
352     
353      DO jj = 1, nlcj
354         DO ji = 1, nlci
355            icof(ji,jj) = idata( mig(ji), mjg(jj) )
356         END DO
357      END DO
358      DO jj = nlcj+1, jpj
359         DO ji = 1, nlci
360            icof(ji,jj) = icof(ji,nlcj)
361         END DO
362      END DO
363      DO jj = 1, jpj
364         DO ji = nlci+1, jpi
365            icof(ji,jj) = icof(nlci,jj)
366         END DO
367      END DO
368
369 9101 FORMAT(1x,a15,2i8)
370 9201 FORMAT(3x,13(i3,12x))
371 9202 FORMAT(i3,41i3)
372
373
374      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
375      ! =================
376      ! define ahm1 and ahm2 at the right grid point position
377      ! (USER: modify ahm1 and ahm2 following your desiderata)
378     
379     
380      ! Decrease ahm to zahmeq m2/s in the tropics
381      ! (from 90   to 20   degrees: ahm = scaled by local metrics
382      !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2
383      !  from  2.5 to  0   degrees: ahm = constant
384      ! symmetric in the south hemisphere)
385
386      zahmeq = aht0
387      zam20s = ahm0*COS( rad * 20. )
388     
389      DO jj = 1, jpj
390         DO ji = 1, jpi
391            IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
392!              leave as set in ldf_dyn_c2d
393            ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
394               ahm2(ji,jj) =  zahmeq
395            ELSE
396               ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
397                  * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
398            ENDIF
399            IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
400!             leave as set in ldf_dyn_c2d
401            ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
402               ahm1(ji,jj) =  zahmeq
403            ELSE
404               ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
405                  * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
406            ENDIF
407         END DO
408      END DO
409
410      ! increase along western boundaries of equatorial strip
411      ! t-point
412      DO jj = 1, jpjm1
413         DO ji = 1, jpim1
414          IF( ABS( gphit(ji,jj) ) < 20. ) THEN
415            zcoft = FLOAT( icof(ji,jj) ) / 100.
416            ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
417          ENDIF
418         END DO
419      END DO
420      ! f-point
421      icof(:,:) = icof(:,:) * tmask(:,:,1)
422      DO jj = 1, jpjm1
423         DO ji = 1, jpim1
424          IF( ABS( gphif(ji,jj) ) < 20. ) THEN
425            zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
426            IF( zmsk == 0. ) THEN
427               zcoff = 1.
428            ELSE
429               zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
430                     / (zmsk * 100.)
431            ENDIF
432            ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
433          ENDIF
434         END DO
435      END DO
436#endif
437     
438      ! Lateral boundary conditions on ( ahm1, ahm2 )
439      !                                ==============
440      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
441      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
442
443      ! Control print
444      IF( lwp .AND. ld_print ) THEN
445         WRITE(numout,*)
446         WRITE(numout,*) 'inildf: 2D ahm1 array'
447         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
448         WRITE(numout,*)
449         WRITE(numout,*) 'inildf: 2D ahm2 array'
450         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
451      ENDIF
452
453      IF( iwrk_not_released(2, 1) )   CALL ctl_stop('ldf_dyn_c2d_orca_R1: failed to release workspace array')
454      !
455   END SUBROUTINE ldf_dyn_c2d_orca_R1
Note: See TracBrowser for help on using the repository browser.