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

Last change on this file since 2690 was 2690, checked in by gm, 12 years ago

dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation

  • 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      END IF
285      !
286   END SUBROUTINE ldf_dyn_c2d_orca
287
288
289   SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print )
290      !!----------------------------------------------------------------------
291      !!                 ***  ROUTINE ldf_dyn_c2d  ***
292      !!
293      !!                   **** W A R N I N G ****
294      !!
295      !!                ORCA R1 configuration
296      !!                 
297      !!                   **** W A R N I N G ****
298      !!                 
299      !! ** Purpose :   initializations of the lateral viscosity for orca R1
300      !!
301      !! ** Method  :   blah blah blah...
302      !!
303      !!----------------------------------------------------------------------
304      USE ldftra_oce, ONLY:   aht0
305      USE wrk_nemo  , ONLY:   iwrk_in_use, iwrk_not_released
306      USE wrk_nemo  , ONLY:   icof => iwrk_2d_1
307      !
308      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
309      !
310      INTEGER ::   ji, jj, jn      ! dummy loop indices
311      INTEGER ::   inum            ! temporary logical unit
312      INTEGER ::   iim, ijm
313      INTEGER ::   ifreq, il1, il2, ij, ii
314      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s
315      CHARACTER (len=15) ::   clexp
316      INTEGER, DIMENSION(jpidta,jpidta) ::   idata
317      !!----------------------------------------------------------------------
318
319      IF( iwrk_in_use(2, 1) ) THEN
320         CALL ctl_stop('ldf_dyn_c2d_orca_R1: requested workspace array is unavailable')   ;   RETURN
321      ENDIF
322
323      IF(lwp) WRITE(numout,*)
324      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
325      IF(lwp) WRITE(numout,*) '~~~~~~  --'
326      IF(lwp) WRITE(numout,*) '        orca_r1 configuration'
327
328#if defined key_antarctic
329#     include "ldfdyn_antarctic.h90"
330#elif defined key_arctic
331#     include "ldfdyn_arctic.h90"
332#else
333      ! Read 2d integer array to specify western boundary increase in the
334      ! ===================== equatorial strip (20N-20S) defined at t-points
335
336      CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
337         &           1, numout, lwp )
338      REWIND inum
339      READ(inum,9101) clexp, iim, ijm
340      READ(inum,'(/)')
341      ifreq = 40
342      il1 = 1
343      DO jn = 1, jpidta/ifreq+1
344         READ(inum,'(/)')
345         il2 = MIN( jpidta, il1+ifreq-1 )
346         READ(inum,9201) ( ii, ji = il1, il2, 5 )
347         READ(inum,'(/)')
348         DO jj = jpjdta, 1, -1
349            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
350         END DO
351         il1 = il1 + ifreq
352      END DO
353     
354      DO jj = 1, nlcj
355         DO ji = 1, nlci
356            icof(ji,jj) = idata( mig(ji), mjg(jj) )
357         END DO
358      END DO
359      DO jj = nlcj+1, jpj
360         DO ji = 1, nlci
361            icof(ji,jj) = icof(ji,nlcj)
362         END DO
363      END DO
364      DO jj = 1, jpj
365         DO ji = nlci+1, jpi
366            icof(ji,jj) = icof(nlci,jj)
367         END DO
368      END DO
369
370 9101 FORMAT(1x,a15,2i8)
371 9201 FORMAT(3x,13(i3,12x))
372 9202 FORMAT(i3,41i3)
373
374
375      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
376      ! =================
377      ! define ahm1 and ahm2 at the right grid point position
378      ! (USER: modify ahm1 and ahm2 following your desiderata)
379     
380     
381      ! Decrease ahm to zahmeq m2/s in the tropics
382      ! (from 90   to 20   degrees: ahm = scaled by local metrics
383      !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2
384      !  from  2.5 to  0   degrees: ahm = constant
385      ! symmetric in the south hemisphere)
386
387      zahmeq = aht0
388      zam20s = ahm0*COS( rad * 20. )
389     
390      DO jj = 1, jpj
391         DO ji = 1, jpi
392            IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
393!              leave as set in ldf_dyn_c2d
394            ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
395               ahm2(ji,jj) =  zahmeq
396            ELSE
397               ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
398                  * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
399            ENDIF
400            IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
401!             leave as set in ldf_dyn_c2d
402            ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
403               ahm1(ji,jj) =  zahmeq
404            ELSE
405               ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
406                  * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
407            ENDIF
408         END DO
409      END DO
410
411      ! increase along western boundaries of equatorial strip
412      ! t-point
413      DO jj = 1, jpjm1
414         DO ji = 1, jpim1
415          IF( ABS( gphit(ji,jj) ) < 20. ) THEN
416            zcoft = FLOAT( icof(ji,jj) ) / 100.
417            ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
418          ENDIF
419         END DO
420      END DO
421      ! f-point
422      icof(:,:) = icof(:,:) * tmask(:,:,1)
423      DO jj = 1, jpjm1
424         DO ji = 1, jpim1
425          IF( ABS( gphif(ji,jj) ) < 20. ) THEN
426            zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
427            IF( zmsk == 0. ) THEN
428               zcoff = 1.
429            ELSE
430               zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
431                     / (zmsk * 100.)
432            ENDIF
433            ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
434          ENDIF
435         END DO
436      END DO
437#endif
438     
439      ! Lateral boundary conditions on ( ahm1, ahm2 )
440      !                                ==============
441      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
442      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
443
444      ! Control print
445      IF( lwp .AND. ld_print ) THEN
446         WRITE(numout,*)
447         WRITE(numout,*) 'inildf: 2D ahm1 array'
448         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
449         WRITE(numout,*)
450         WRITE(numout,*) 'inildf: 2D ahm2 array'
451         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
452      ENDIF
453
454      IF( iwrk_not_released(2, 1) )   CALL ctl_stop('ldf_dyn_c2d_orca_R1: failed to release workspace array')
455      !
456   END SUBROUTINE ldf_dyn_c2d_orca_R1
Note: See TracBrowser for help on using the repository browser.