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

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 12555

Last change on this file since 12555 was 12555, checked in by charris, 4 years ago

Changes from GO6 package branch (GMED ticket 450):

svn merge -r 11035:11101 svn+ssh://charris@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package

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.