source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 18 months ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 20.6 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      IF(lwp .AND. lflush) CALL flush(numout)
44
45      ! harmonic operator (ahm1, ahm2) : ( T- and F- points) (used for laplacian operators
46      ! ===============================                       whatever its orientation is)
47      IF( ln_dynldf_lap ) THEN
48         ! define ahm1 and ahm2 at the right grid point position
49         ! (USER: modify ahm1 and ahm2 following your desiderata)
50
51         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
52         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
53
54         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
55         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
56
57         za00 = ahm0 / zd_max
58         DO jj = 1, jpj
59            DO ji = 1, jpi
60               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
61               zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
62               ahm1(ji,jj) = za00 * zetmax
63               ahm2(ji,jj) = za00 * zefmax
64            END DO
65         END DO
66
67         IF( ln_dynldf_iso ) THEN
68            IF(lwp) WRITE(numout,*) '              Caution, as implemented now, the isopycnal part of momentum'
69            IF(lwp) WRITE(numout,*) '                 mixing use aht0 as eddy viscosity coefficient. Thus, it is'
70            IF(lwp) WRITE(numout,*) '                 uniform and you must be sure that your ahm is greater than'
71            IF(lwp) WRITE(numout,*) '                 aht0 everywhere in the model domain.'
72         ENDIF
73
74         ! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2)
75         ! ==============================================
76         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   CALL ldf_dyn_c2d_orca( ld_print )
77         IF( cp_cfg == "orca" .AND.   jp_cfg == 1)                       CALL ldf_dyn_c2d_orca_R1( 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         IF(lwp .AND. lflush) CALL flush(numout)
89      ENDIF
90
91      ! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator
92      ! =================================                      whatever its orientation is)
93      IF( ln_dynldf_bilap ) THEN
94         ! (USER: modify ahm3 and ahm4 following your desiderata)
95         ! Here: ahm is proportional to the cube of the maximum of the gridspacing
96         !       in the to horizontal direction
97
98         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
99         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
100
101         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
102         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
103
104         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
105         DO jj = 1, jpj
106            DO ji = 1, jpi
107               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
108               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
109               ahm3(ji,jj) = za00 * zeumax * zeumax * zeumax
110               ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax
111            END DO
112         END DO
113
114         ! Control print
115         IF( lwp .AND. ld_print ) THEN
116            WRITE(numout,*)
117            WRITE(numout,*) 'inildf: ahm3 array'
118            CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
119            WRITE(numout,*)
120            WRITE(numout,*) 'inildf: ahm4 array'
121            CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
122         ENDIF
123         IF(lwp .AND. lflush) CALL flush(numout)
124      ENDIF
125      !
126   END SUBROUTINE ldf_dyn_c2d
127
128
129   SUBROUTINE ldf_dyn_c2d_orca( ld_print )
130      !!----------------------------------------------------------------------
131      !!                 ***  ROUTINE ldf_dyn_c2d  ***
132      !!
133      !!                   **** W A R N I N G ****
134      !!
135      !!                ORCA R2 and R4 configurations
136      !!                 
137      !!                   **** W A R N I N G ****
138      !!                 
139      !! ** Purpose :   initializations of the lateral viscosity for orca R2
140      !!
141      !! ** Method  :   blah blah blah...
142      !!
143      !!----------------------------------------------------------------------
144      USE ldftra_oce, ONLY:   aht0
145      USE iom
146      !
147      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
148      !
149      INTEGER  ::   ji, jj, jn   ! dummy loop indices
150      INTEGER  ::   inum, iim, ijm            ! local integers
151      INTEGER  ::   ifreq, il1, il2, ij, ii
152      INTEGER  ::   ijpt0,ijpt1, ierror
153      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
154      CHARACTER (len=15) ::   clexp
155      INTEGER,     POINTER, DIMENSION(:,:)  :: icof
156      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztemp2d  ! temporary array to read ahmcoef file
157      !!----------------------------------------------------------------------
158      !                               
159      CALL wrk_alloc( jpi   , jpj   , icof  )
160      !
161      IF(lwp) WRITE(numout,*)
162      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
163      IF(lwp) WRITE(numout,*) '~~~~~~  --'
164      IF(lwp) WRITE(numout,*) '        orca ocean configuration'
165      IF(lwp .AND. lflush) CALL flush(numout)
166
167      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
168!
169! 1.2 Modify ahm
170! --------------
171         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
172         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
173         IF(lwp)WRITE(numout,*) '         north boundary increase'
174
175         ahm1(:,:) = ahm0
176         ahm2(:,:) = ahm0
177
178         ijpt0=max(1,min(49 -njmpp+1,jpj))
179         ijpt1=max(0,min(49-njmpp+1,jpj-1))
180         DO jj=ijpt0,ijpt1
181            ahm2(:,jj)=ahm0*2.
182            ahm1(:,jj)=ahm0*2.
183         END DO
184         ijpt0=max(1,min(48 -njmpp+1,jpj))
185         ijpt1=max(0,min(48-njmpp+1,jpj-1))
186         DO jj=ijpt0,ijpt1
187            ahm2(:,jj)=ahm0*1.9
188            ahm1(:,jj)=ahm0*1.75
189         END DO
190         ijpt0=max(1,min(47 -njmpp+1,jpj))
191         ijpt1=max(0,min(47-njmpp+1,jpj-1))
192         DO jj=ijpt0,ijpt1
193            ahm2(:,jj)=ahm0*1.5
194            ahm1(:,jj)=ahm0*1.25
195         END DO
196         ijpt0=max(1,min(46 -njmpp+1,jpj))
197         ijpt1=max(0,min(46-njmpp+1,jpj-1))
198         DO jj=ijpt0,ijpt1
199            ahm2(:,jj)=ahm0*1.1
200         END DO
201
202      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
203! 1.2 Modify ahm
204! --------------
205         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
206         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
207         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
208
209
210         ahm1(:,:) = ahm0
211         ahm2(:,:) = ahm0
212
213         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
214         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
215         DO jj=ijpt0,ijpt1
216            ahm2(:,jj)=ahm0*2.
217            ahm1(:,jj)=ahm0*2.
218         END DO
219         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
220         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
221         DO jj=ijpt0,ijpt1
222            ahm2(:,jj)=ahm0*1.9
223            ahm1(:,jj)=ahm0*1.75
224         END DO
225         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
226         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
227         DO jj=ijpt0,ijpt1
228            ahm2(:,jj)=ahm0*1.5
229            ahm1(:,jj)=ahm0*1.25
230         END DO
231         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
232         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
233         DO jj=ijpt0,ijpt1
234            ahm2(:,jj)=ahm0*1.1
235         END DO
236      ELSE
237         ! Read 2d integer array to specify western boundary increase in the
238         ! ===================== equatorial strip (20N-20S) defined at t-points
239         !
240         ALLOCATE( ztemp2d(jpi,jpj) )
241         ztemp2d(:,:) = 0.
242         CALL iom_open ( 'ahmcoef.nc', inum )
243         CALL iom_get  ( inum, jpdom_data, 'icof', ztemp2d)
244         icof(:,:)  = NINT(ztemp2d(:,:))
245         CALL iom_close( inum )
246         DEALLOCATE(ztemp2d)
247
248         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
249         ! =================
250         ! define ahm1 and ahm2 at the right grid point position
251         ! (USER: modify ahm1 and ahm2 following your desiderata)
252
253
254         ! Decrease ahm to zahmeq m2/s in the tropics
255         ! (from 90 to 20 degre: ahm = constant
256         ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
257         ! from  2.5 to  0 degre: ahm = constant
258         ! symmetric in the south hemisphere)
259
260         zahmeq = aht0
261
262         DO jj = 1, jpj
263            DO ji = 1, jpi
264               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
265                  ahm2(ji,jj) =  ahm0
266               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
267                  ahm2(ji,jj) =  zahmeq
268               ELSE
269                  ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
270                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
271               ENDIF
272               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
273                  ahm1(ji,jj) =  ahm0
274               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
275                  ahm1(ji,jj) =  zahmeq
276               ELSE
277                  ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
278                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
279               ENDIF
280            END DO
281         END DO
282
283         ! increase along western boundaries of equatorial strip
284         ! t-point
285         DO jj = 1, jpjm1
286            DO ji = 1, jpim1
287               zcoft = FLOAT( icof(ji,jj) ) / 100.
288               ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
289            END DO
290         END DO
291         ! f-point
292         icof(:,:) = icof(:,:) * tmask(:,:,1)
293         DO jj = 1, jpjm1
294            DO ji = 1, jpim1   ! NO vector opt.
295               zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
296               IF( zmsk == 0. ) THEN
297                  zcoff = 1.
298               ELSE
299                  zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
300                     / (zmsk * 100.)
301               ENDIF
302               ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
303            END DO
304         END DO
305      ENDIF
306     
307      ! Lateral boundary conditions on ( ahm1, ahm2 )
308      !                                ==============
309      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
310      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
311
312      ! Control print
313      IF( lwp .AND. ld_print ) THEN
314         WRITE(numout,*)
315         WRITE(numout,*) 'inildf: 2D ahm1 array'
316         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
317         WRITE(numout,*)
318         WRITE(numout,*) 'inildf: 2D ahm2 array'
319         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
320      ENDIF
321      !
322      IF(lwp .AND. lflush) CALL flush(numout)
323      !
324      CALL wrk_dealloc( jpi   , jpj   , icof  )
325      !
326   END SUBROUTINE ldf_dyn_c2d_orca
327
328
329   SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print )
330      !!----------------------------------------------------------------------
331      !!                 ***  ROUTINE ldf_dyn_c2d  ***
332      !!
333      !!                   **** W A R N I N G ****
334      !!
335      !!                ORCA R1 configuration
336      !!                 
337      !!                   **** W A R N I N G ****
338      !!                 
339      !! ** Purpose :   initializations of the lateral viscosity for orca R1
340      !!
341      !! ** Method  :   blah blah blah...
342      !!
343      !!----------------------------------------------------------------------
344      USE ldftra_oce, ONLY:   aht0
345      USE iom
346      !
347      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
348      !
349      INTEGER ::   ji, jj, jn      ! dummy loop indices
350      INTEGER ::   inum            ! temporary logical unit
351      INTEGER ::   iim, ijm
352      INTEGER ::   ifreq, il1, il2, ij, ii
353      INTEGER ::   ijpt0,ijpt1, ierror
354      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s
355      CHARACTER (len=15) ::   clexp
356      INTEGER,     POINTER, DIMENSION(:,:)  :: icof
357      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztemp2d  ! temporary array to read ahmcoef file
358      !!----------------------------------------------------------------------
359      !                               
360      CALL wrk_alloc( jpi   , jpj   , icof  )
361      !                               
362      IF(lwp) WRITE(numout,*)
363      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
364      IF(lwp) WRITE(numout,*) '~~~~~~  --'
365      IF(lwp) WRITE(numout,*) '        orca_r1 configuration'
366      IF(lwp .AND. lflush) CALL flush(numout)
367
368      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
369!
370! 1.2 Modify ahm
371! --------------
372         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
373         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
374         IF(lwp)WRITE(numout,*) '         north boundary increase'
375
376         ahm1(:,:) = ahm0
377         ahm2(:,:) = ahm0
378
379         ijpt0=max(1,min(49 -njmpp+1,jpj))
380         ijpt1=max(0,min(49-njmpp+1,jpj-1))
381         DO jj=ijpt0,ijpt1
382            ahm2(:,jj)=ahm0*2.
383            ahm1(:,jj)=ahm0*2.
384         END DO
385         ijpt0=max(1,min(48 -njmpp+1,jpj))
386         ijpt1=max(0,min(48-njmpp+1,jpj-1))
387         DO jj=ijpt0,ijpt1
388            ahm2(:,jj)=ahm0*1.9
389            ahm1(:,jj)=ahm0*1.75
390         END DO
391         ijpt0=max(1,min(47 -njmpp+1,jpj))
392         ijpt1=max(0,min(47-njmpp+1,jpj-1))
393         DO jj=ijpt0,ijpt1
394            ahm2(:,jj)=ahm0*1.5
395            ahm1(:,jj)=ahm0*1.25
396         END DO
397         ijpt0=max(1,min(46 -njmpp+1,jpj))
398         ijpt1=max(0,min(46-njmpp+1,jpj-1))
399         DO jj=ijpt0,ijpt1
400            ahm2(:,jj)=ahm0*1.1
401         END DO
402
403      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
404! 1.2 Modify ahm
405! --------------
406         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
407         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
408         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
409
410
411         ahm1(:,:) = ahm0
412         ahm2(:,:) = ahm0
413
414         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
415         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
416         DO jj=ijpt0,ijpt1
417            ahm2(:,jj)=ahm0*2.
418            ahm1(:,jj)=ahm0*2.
419         END DO
420         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
421         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
422         DO jj=ijpt0,ijpt1
423            ahm2(:,jj)=ahm0*1.9
424            ahm1(:,jj)=ahm0*1.75
425         END DO
426         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
427         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
428         DO jj=ijpt0,ijpt1
429            ahm2(:,jj)=ahm0*1.5
430            ahm1(:,jj)=ahm0*1.25
431         END DO
432         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
433         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
434         DO jj=ijpt0,ijpt1
435            ahm2(:,jj)=ahm0*1.1
436         END DO
437      ELSE
438         
439         ! Read 2d integer array to specify western boundary increase in the
440         ! ===================== equatorial strip (20N-20S) defined at t-points
441         ALLOCATE( ztemp2d(jpi,jpj) )
442         ztemp2d(:,:) = 0.
443         CALL iom_open ( 'ahmcoef.nc', inum )
444         CALL iom_get  ( inum, jpdom_data, 'icof', ztemp2d)
445         icof(:,:)  = NINT(ztemp2d(:,:))
446         CALL iom_close( inum )
447         DEALLOCATE(ztemp2d)
448
449         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
450         ! =================
451         ! define ahm1 and ahm2 at the right grid point position
452         ! (USER: modify ahm1 and ahm2 following your desiderata)
453
454
455         ! Decrease ahm to zahmeq m2/s in the tropics
456         ! (from 90   to 20   degrees: ahm = scaled by local metrics
457         !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2
458         !  from  2.5 to  0   degrees: ahm = constant
459         ! symmetric in the south hemisphere)
460
461         zahmeq = aht0
462         zam20s = ahm0*COS( rad * 20. )
463
464         DO jj = 1, jpj
465            DO ji = 1, jpi
466               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
467                  !              leave as set in ldf_dyn_c2d
468               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
469                  ahm2(ji,jj) =  zahmeq
470               ELSE
471                  ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
472                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
473               ENDIF
474               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
475                  !             leave as set in ldf_dyn_c2d
476               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
477                  ahm1(ji,jj) =  zahmeq
478               ELSE
479                  ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
480                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
481               ENDIF
482            END DO
483         END DO
484
485         ! increase along western boundaries of equatorial strip
486         ! t-point
487         DO jj = 1, jpjm1
488            DO ji = 1, jpim1
489               IF( ABS( gphit(ji,jj) ) < 20. ) THEN
490                  zcoft = FLOAT( icof(ji,jj) ) / 100.
491                  ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
492               ENDIF
493            END DO
494         END DO
495         ! f-point
496         icof(:,:) = icof(:,:) * tmask(:,:,1)
497         DO jj = 1, jpjm1
498            DO ji = 1, jpim1
499               IF( ABS( gphif(ji,jj) ) < 20. ) THEN
500                  zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
501                  IF( zmsk == 0. ) THEN
502                     zcoff = 1.
503                  ELSE
504                     zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
505                        / (zmsk * 100.)
506                  ENDIF
507                  ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
508               ENDIF
509            END DO
510         END DO
511      ENDIF
512     
513      ! Lateral boundary conditions on ( ahm1, ahm2 )
514      !                                ==============
515      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
516      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
517
518      ! Control print
519      IF( lwp .AND. ld_print ) THEN
520         WRITE(numout,*)
521         WRITE(numout,*) 'inildf: 2D ahm1 array'
522         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
523         WRITE(numout,*)
524         WRITE(numout,*) 'inildf: 2D ahm2 array'
525         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
526      ENDIF
527      !
528      IF(lwp .AND. lflush) CALL flush(numout)
529      !
530      CALL wrk_dealloc( jpi   , jpj   , icof  )
531      !
532   END SUBROUTINE ldf_dyn_c2d_orca_R1
Note: See TracBrowser for help on using the repository browser.