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

source: trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 5040

Last change on this file since 5040 was 5040, checked in by smasson, 10 years ago

supress wrk_alloc( jpidta, jpjdta, see ticket #1440
  • Property svn:keywords set to Id
File size: 22.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      !
143      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
144      !
145      INTEGER  ::   ji, jj, jn   ! dummy loop indices
146      INTEGER  ::   inum, iim, ijm            ! local integers
147      INTEGER  ::   ifreq, il1, il2, ij, ii
148      INTEGER  ::   ijpt0,ijpt1, ierror
149      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
150      CHARACTER (len=15) ::   clexp
151      INTEGER,     POINTER, DIMENSION(:,:)  :: icof
152      INTEGER, ALLOCATABLE, DIMENSION(:,:)  :: idata
153      !!----------------------------------------------------------------------
154      !                               
155      CALL wrk_alloc( jpi   , jpj   , icof  )
156      !
157      IF(lwp) WRITE(numout,*)
158      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
159      IF(lwp) WRITE(numout,*) '~~~~~~  --'
160      IF(lwp) WRITE(numout,*) '        orca ocean configuration'
161
162      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
163!
164! 1.2 Modify ahm
165! --------------
166         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
167         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
168         IF(lwp)WRITE(numout,*) '         north boundary increase'
169
170         ahm1(:,:) = ahm0
171         ahm2(:,:) = ahm0
172
173         ijpt0=max(1,min(49 -njmpp+1,jpj))
174         ijpt1=max(0,min(49-njmpp+1,jpj-1))
175         DO jj=ijpt0,ijpt1
176            ahm2(:,jj)=ahm0*2.
177            ahm1(:,jj)=ahm0*2.
178         END DO
179         ijpt0=max(1,min(48 -njmpp+1,jpj))
180         ijpt1=max(0,min(48-njmpp+1,jpj-1))
181         DO jj=ijpt0,ijpt1
182            ahm2(:,jj)=ahm0*1.9
183            ahm1(:,jj)=ahm0*1.75
184         END DO
185         ijpt0=max(1,min(47 -njmpp+1,jpj))
186         ijpt1=max(0,min(47-njmpp+1,jpj-1))
187         DO jj=ijpt0,ijpt1
188            ahm2(:,jj)=ahm0*1.5
189            ahm1(:,jj)=ahm0*1.25
190         END DO
191         ijpt0=max(1,min(46 -njmpp+1,jpj))
192         ijpt1=max(0,min(46-njmpp+1,jpj-1))
193         DO jj=ijpt0,ijpt1
194            ahm2(:,jj)=ahm0*1.1
195         END DO
196
197      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
198! 1.2 Modify ahm
199! --------------
200         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
201         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
202         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
203
204
205         ahm1(:,:) = ahm0
206         ahm2(:,:) = ahm0
207
208         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
209         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
210         DO jj=ijpt0,ijpt1
211            ahm2(:,jj)=ahm0*2.
212            ahm1(:,jj)=ahm0*2.
213         END DO
214         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
215         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
216         DO jj=ijpt0,ijpt1
217            ahm2(:,jj)=ahm0*1.9
218            ahm1(:,jj)=ahm0*1.75
219         END DO
220         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
221         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
222         DO jj=ijpt0,ijpt1
223            ahm2(:,jj)=ahm0*1.5
224            ahm1(:,jj)=ahm0*1.25
225         END DO
226         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
227         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
228         DO jj=ijpt0,ijpt1
229            ahm2(:,jj)=ahm0*1.1
230         END DO
231      ELSE
232         ! Read 2d integer array to specify western boundary increase in the
233         ! ===================== equatorial strip (20N-20S) defined at t-points
234         
235         ALLOCATE( idata(jpidta,jpjdta), STAT=ierror )
236         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_c2d_orca: unable to allocate idata array' )
237         !
238         CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
239         READ(inum,9101) clexp, iim, ijm
240         READ(inum,'(/)')
241         ifreq = 40
242         il1 = 1
243         DO jn = 1, jpidta/ifreq+1
244            READ(inum,'(/)')
245            il2 = MIN( jpidta, il1+ifreq-1 )
246            READ(inum,9201) ( ii, ji = il1, il2, 5 )
247            READ(inum,'(/)')
248            DO jj = jpjdta, 1, -1
249               READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
250            END DO
251            il1 = il1 + ifreq
252         END DO
253
254         DO jj = 1, nlcj
255            DO ji = 1, nlci
256               icof(ji,jj) = idata( mig(ji), mjg(jj) )
257            END DO
258         END DO
259         DO jj = nlcj+1, jpj
260            DO ji = 1, nlci
261               icof(ji,jj) = icof(ji,nlcj)
262            END DO
263         END DO
264         DO jj = 1, jpj
265            DO ji = nlci+1, jpi
266               icof(ji,jj) = icof(nlci,jj)
267            END DO
268         END DO
269
2709101     FORMAT(1x,a15,2i8)
2719201     FORMAT(3x,13(i3,12x))
2729202     FORMAT(i3,41i3)
273         
274         DEALLOCATE(idata)
275
276         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
277         ! =================
278         ! define ahm1 and ahm2 at the right grid point position
279         ! (USER: modify ahm1 and ahm2 following your desiderata)
280
281
282         ! Decrease ahm to zahmeq m2/s in the tropics
283         ! (from 90 to 20 degre: ahm = constant
284         ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
285         ! from  2.5 to  0 degre: ahm = constant
286         ! symmetric in the south hemisphere)
287
288         zahmeq = aht0
289
290         DO jj = 1, jpj
291            DO ji = 1, jpi
292               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
293                  ahm2(ji,jj) =  ahm0
294               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
295                  ahm2(ji,jj) =  zahmeq
296               ELSE
297                  ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
298                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
299               ENDIF
300               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
301                  ahm1(ji,jj) =  ahm0
302               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
303                  ahm1(ji,jj) =  zahmeq
304               ELSE
305                  ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
306                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
307               ENDIF
308            END DO
309         END DO
310
311         ! increase along western boundaries of equatorial strip
312         ! t-point
313         DO jj = 1, jpjm1
314            DO ji = 1, jpim1
315               zcoft = FLOAT( icof(ji,jj) ) / 100.
316               ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
317            END DO
318         END DO
319         ! f-point
320         icof(:,:) = icof(:,:) * tmask(:,:,1)
321         DO jj = 1, jpjm1
322            DO ji = 1, jpim1   ! NO vector opt.
323               zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
324               IF( zmsk == 0. ) THEN
325                  zcoff = 1.
326               ELSE
327                  zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
328                     / (zmsk * 100.)
329               ENDIF
330               ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
331            END DO
332         END DO
333      ENDIF
334     
335      ! Lateral boundary conditions on ( ahm1, ahm2 )
336      !                                ==============
337      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
338      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
339
340      ! Control print
341      IF( lwp .AND. ld_print ) THEN
342         WRITE(numout,*)
343         WRITE(numout,*) 'inildf: 2D ahm1 array'
344         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
345         WRITE(numout,*)
346         WRITE(numout,*) 'inildf: 2D ahm2 array'
347         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
348      ENDIF
349      !
350      CALL wrk_dealloc( jpi   , jpj   , icof  )
351      !
352   END SUBROUTINE ldf_dyn_c2d_orca
353
354
355   SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print )
356      !!----------------------------------------------------------------------
357      !!                 ***  ROUTINE ldf_dyn_c2d  ***
358      !!
359      !!                   **** W A R N I N G ****
360      !!
361      !!                ORCA R1 configuration
362      !!                 
363      !!                   **** W A R N I N G ****
364      !!                 
365      !! ** Purpose :   initializations of the lateral viscosity for orca R1
366      !!
367      !! ** Method  :   blah blah blah...
368      !!
369      !!----------------------------------------------------------------------
370      USE ldftra_oce, ONLY:   aht0
371      !
372      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
373      !
374      INTEGER ::   ji, jj, jn      ! dummy loop indices
375      INTEGER ::   inum            ! temporary logical unit
376      INTEGER ::   iim, ijm
377      INTEGER ::   ifreq, il1, il2, ij, ii
378      INTEGER ::   ijpt0,ijpt1, ierror
379      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s
380      CHARACTER (len=15) ::   clexp
381      INTEGER,     POINTER, DIMENSION(:,:)  :: icof
382      INTEGER, ALLOCATABLE, DIMENSION(:,:)  :: idata
383      !!----------------------------------------------------------------------
384      !                               
385      CALL wrk_alloc( jpi   , jpj   , icof  )
386      !                               
387
388      IF(lwp) WRITE(numout,*)
389      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
390      IF(lwp) WRITE(numout,*) '~~~~~~  --'
391      IF(lwp) WRITE(numout,*) '        orca_r1 configuration'
392
393      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
394!
395! 1.2 Modify ahm
396! --------------
397         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
398         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
399         IF(lwp)WRITE(numout,*) '         north boundary increase'
400
401         ahm1(:,:) = ahm0
402         ahm2(:,:) = ahm0
403
404         ijpt0=max(1,min(49 -njmpp+1,jpj))
405         ijpt1=max(0,min(49-njmpp+1,jpj-1))
406         DO jj=ijpt0,ijpt1
407            ahm2(:,jj)=ahm0*2.
408            ahm1(:,jj)=ahm0*2.
409         END DO
410         ijpt0=max(1,min(48 -njmpp+1,jpj))
411         ijpt1=max(0,min(48-njmpp+1,jpj-1))
412         DO jj=ijpt0,ijpt1
413            ahm2(:,jj)=ahm0*1.9
414            ahm1(:,jj)=ahm0*1.75
415         END DO
416         ijpt0=max(1,min(47 -njmpp+1,jpj))
417         ijpt1=max(0,min(47-njmpp+1,jpj-1))
418         DO jj=ijpt0,ijpt1
419            ahm2(:,jj)=ahm0*1.5
420            ahm1(:,jj)=ahm0*1.25
421         END DO
422         ijpt0=max(1,min(46 -njmpp+1,jpj))
423         ijpt1=max(0,min(46-njmpp+1,jpj-1))
424         DO jj=ijpt0,ijpt1
425            ahm2(:,jj)=ahm0*1.1
426         END DO
427
428      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
429! 1.2 Modify ahm
430! --------------
431         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
432         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
433         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
434
435
436         ahm1(:,:) = ahm0
437         ahm2(:,:) = ahm0
438
439         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
440         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
441         DO jj=ijpt0,ijpt1
442            ahm2(:,jj)=ahm0*2.
443            ahm1(:,jj)=ahm0*2.
444         END DO
445         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
446         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
447         DO jj=ijpt0,ijpt1
448            ahm2(:,jj)=ahm0*1.9
449            ahm1(:,jj)=ahm0*1.75
450         END DO
451         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
452         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
453         DO jj=ijpt0,ijpt1
454            ahm2(:,jj)=ahm0*1.5
455            ahm1(:,jj)=ahm0*1.25
456         END DO
457         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
458         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
459         DO jj=ijpt0,ijpt1
460            ahm2(:,jj)=ahm0*1.1
461         END DO
462      ELSE
463         
464         ! Read 2d integer array to specify western boundary increase in the
465         ! ===================== equatorial strip (20N-20S) defined at t-points
466         
467         ALLOCATE( idata(jpidta,jpjdta), STAT=ierror )
468         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_c2d_orca_R1: unable to allocate idata array' )
469         !
470         CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
471            &           1, numout, lwp )
472         REWIND inum
473         READ(inum,9101) clexp, iim, ijm
474         READ(inum,'(/)')
475         ifreq = 40
476         il1 = 1
477         DO jn = 1, jpidta/ifreq+1
478            READ(inum,'(/)')
479            il2 = MIN( jpidta, il1+ifreq-1 )
480            READ(inum,9201) ( ii, ji = il1, il2, 5 )
481            READ(inum,'(/)')
482            DO jj = jpjdta, 1, -1
483               READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
484            END DO
485            il1 = il1 + ifreq
486         END DO
487
488         DO jj = 1, nlcj
489            DO ji = 1, nlci
490               icof(ji,jj) = idata( mig(ji), mjg(jj) )
491            END DO
492         END DO
493         DO jj = nlcj+1, jpj
494            DO ji = 1, nlci
495               icof(ji,jj) = icof(ji,nlcj)
496            END DO
497         END DO
498         DO jj = 1, jpj
499            DO ji = nlci+1, jpi
500               icof(ji,jj) = icof(nlci,jj)
501            END DO
502         END DO
503
5049101     FORMAT(1x,a15,2i8)
5059201     FORMAT(3x,13(i3,12x))
5069202     FORMAT(i3,41i3)
507         
508         DEALLOCATE(idata)
509
510         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
511         ! =================
512         ! define ahm1 and ahm2 at the right grid point position
513         ! (USER: modify ahm1 and ahm2 following your desiderata)
514
515
516         ! Decrease ahm to zahmeq m2/s in the tropics
517         ! (from 90   to 20   degrees: ahm = scaled by local metrics
518         !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2
519         !  from  2.5 to  0   degrees: ahm = constant
520         ! symmetric in the south hemisphere)
521
522         zahmeq = aht0
523         zam20s = ahm0*COS( rad * 20. )
524
525         DO jj = 1, jpj
526            DO ji = 1, jpi
527               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
528                  !              leave as set in ldf_dyn_c2d
529               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
530                  ahm2(ji,jj) =  zahmeq
531               ELSE
532                  ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
533                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
534               ENDIF
535               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
536                  !             leave as set in ldf_dyn_c2d
537               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
538                  ahm1(ji,jj) =  zahmeq
539               ELSE
540                  ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
541                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
542               ENDIF
543            END DO
544         END DO
545
546         ! increase along western boundaries of equatorial strip
547         ! t-point
548         DO jj = 1, jpjm1
549            DO ji = 1, jpim1
550               IF( ABS( gphit(ji,jj) ) < 20. ) THEN
551                  zcoft = FLOAT( icof(ji,jj) ) / 100.
552                  ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
553               ENDIF
554            END DO
555         END DO
556         ! f-point
557         icof(:,:) = icof(:,:) * tmask(:,:,1)
558         DO jj = 1, jpjm1
559            DO ji = 1, jpim1
560               IF( ABS( gphif(ji,jj) ) < 20. ) THEN
561                  zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
562                  IF( zmsk == 0. ) THEN
563                     zcoff = 1.
564                  ELSE
565                     zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
566                        / (zmsk * 100.)
567                  ENDIF
568                  ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
569               ENDIF
570            END DO
571         END DO
572      ENDIF
573     
574      ! Lateral boundary conditions on ( ahm1, ahm2 )
575      !                                ==============
576      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
577      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
578
579      ! Control print
580      IF( lwp .AND. ld_print ) THEN
581         WRITE(numout,*)
582         WRITE(numout,*) 'inildf: 2D ahm1 array'
583         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
584         WRITE(numout,*)
585         WRITE(numout,*) 'inildf: 2D ahm2 array'
586         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
587      ENDIF
588      !
589      CALL wrk_dealloc( jpi   , jpj   , icof  )
590      !
591   END SUBROUTINE ldf_dyn_c2d_orca_R1
Note: See TracBrowser for help on using the repository browser.