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_c3d.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_c3d.h90 @ 2620

Last change on this file since 2620 was 2590, checked in by trackstand2, 13 years ago

Merge branch 'dynamic_memory' into master-svn-dyn
  • Property svn:keywords set to Id
File size: 19.0 KB
Line 
1   !!----------------------------------------------------------------------
2   !!                        ***  ldfdyn_c3d.h90  ***
3   !!----------------------------------------------------------------------
4
5   !!----------------------------------------------------------------------
6   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
7   !! $Id$
8   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   'key_dynldf_c3d'             3D lateral eddy viscosity coefficients
13   !!----------------------------------------------------------------------
14
15   SUBROUTINE ldf_dyn_c3d( ld_print )
16      !!----------------------------------------------------------------------
17      !!                  ***  ROUTINE ldf_dyn_c3d  ***
18      !!                   
19      !! ** Purpose :   initializations of the horizontal ocean physics
20      !!
21      !! ** Method  :   3D eddy viscosity coef. ( longitude, latitude, depth )
22      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
23      !!                              ahm2, ahm4 never used
24      !!       bilaplacian operator : ahm1, ahm2 never used
25      !!                           :  ahm3, ahm4 defined at U- and V-points
26      !!       ??? explanation of the default is missing
27      !!----------------------------------------------------------------------
28      USE ldftra_oce, ONLY : aht0
29      USE wrk_nemo, ONLY: wrk_use, wrk_release
30      USE wrk_nemo, ONLY: zcoef => wrk_1d_2
31      !!
32      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
33      !!
34      INTEGER  ::   ji, jj, jk      ! dummy loop indices
35      REAL(wp) ::   &
36         zr = 0.2 ,   &  ! maximum of the reduction factor at the bottom ocean
37         !               ! ( 0 < zr < 1 )
38         zh = 500.,   &  ! depth of at which start the reduction ( > dept(1) )
39         zd_max   ,   &  ! maximum grid spacing over the global domain
40         za00, zc, zd    ! temporary scalars
41      REAL(wp) ::        &
42         zetmax, zefmax, &
43         zeumax, zevmax   
44      !!----------------------------------------------------------------------
45
46      IF(.not. wrk_use(1,2))THEN
47         CALL ctl_stop('ldf_dyn_c3d: ERROR: requested workspace array unavailable.')
48         RETURN
49      END IF
50
51      IF(lwp) WRITE(numout,*)
52      IF(lwp) WRITE(numout,*) 'ldf_dyn_c3d : 3D lateral eddy viscosity coefficient'
53      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
54
55     
56      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operators
57      ! =================                       whatever its orientation is)
58      IF( ln_dynldf_lap ) THEN
59         ! define ahm1 and ahm2 at the right grid point position
60         ! (USER: modify ahm1 and ahm2 following your desiderata)
61
62         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
63         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
64
65         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
66         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
67
68         za00 = ahm0 / zd_max
69
70         IF( ln_dynldf_iso ) THEN
71            IF(lwp) WRITE(numout,*) '              Caution, as implemented now, the isopycnal part of momentum'
72            IF(lwp) WRITE(numout,*) '                 mixing use aht0 as eddy viscosity coefficient. Thus, it is'
73            IF(lwp) WRITE(numout,*) '                 uniform and you must be sure that your ahm is greater than'
74            IF(lwp) WRITE(numout,*) '                 aht0 everywhere in the model domain.'
75         ENDIF
76
77         CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm1 )   ! vertical profile
78         CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 )   ! vertical profile
79         DO jk = 1,jpk
80            DO jj = 1, jpj
81               DO ji = 1, jpi
82                  zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
83                  zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
84                  ahm1(ji,jj,jk) = za00 * zetmax * ahm1(ji,jj,jk)
85                  ahm2(ji,jj,jk) = za00 * zefmax * ahm2(ji,jj,jk)
86               END DO
87            END DO
88         END DO
89
90
91         ! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2)
92         ! ==============================================
93         IF( cp_cfg == "orca" .AND. ( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) ) THEN
94            IF(lwp) WRITE(numout,*)
95            IF(lwp) WRITE(numout,*) '              ORCA R1, R2 or R4: overwrite the previous definition of ahm'
96            IF(lwp) WRITE(numout,*) '              ================='
97            CALL ldf_dyn_c3d_orca( ld_print )
98         ENDIF
99
100      ENDIF
101     
102      ! Control print
103      IF(lwp .AND. ld_print ) THEN
104         WRITE(numout,*)
105         WRITE(numout,*) '         3D ahm1 array (k=1)'
106         CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
107         WRITE(numout,*)
108         WRITE(numout,*) '         3D ahm2 array (k=1)'
109         CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
110      ENDIF
111
112
113      ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
114      ! ================================  whatever its orientation is)
115      ! (USER: modify ahm3 and ahm4 following your desiderata)
116      ! Here: ahm is proportional to the cube of the maximum of the gridspacing
117      !       in the to horizontal direction
118
119      IF( ln_dynldf_bilap ) THEN
120
121         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
122         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
123
124         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
125         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
126
127         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
128         DO jj = 1, jpj
129            DO ji = 1, jpi
130               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
131               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
132               ahm3(ji,jj,1) = za00 * zeumax * zeumax * zeumax
133               ahm4(ji,jj,1) = za00 * zevmax * zevmax * zevmax
134            END DO
135         END DO
136
137         zh = MAX( zh, fsdept(1,1,1) )   ! at least the first reach ahm0
138         IF( ln_zco ) THEN               ! z-coordinate, same profile everywhere
139            IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
140            DO jk = 1, jpk
141               IF( fsdept(1,1,jk) <= zh ) THEN
142                  zcoef(jk) = 1.e0
143               ELSE
144                  zcoef(jk) = 1.e0 + ( zr - 1.e0 )   &
145                     &               * (  1. - EXP( ( fsdept(1,1,jk   ) - zh ) / zh )  )   &
146                     &               / (  1. - EXP( ( fsdept(1,1,jpkm1) - zh ) / zh )  )
147               ENDIF
148               ahm3(:,:,jk) = ahm3(:,:,1) * zcoef(jk)
149               ahm4(:,:,jk) = ahm4(:,:,1) * zcoef(jk)
150               IF(lwp) WRITE(numout,'(34x,E7.2,8x,i3)') zcoef(jk) * ahm0, jk
151            END DO
152         ELSE                            ! partial steps or s-ccordinate
153            zc = MAXVAL( fsdept(:,:,jpkm1) )
154            IF( lk_mpp )   CALL mpp_max( zc )   ! max over the global domain
155
156            zc = 1. / (  1. - EXP( ( zc - zh ) / zh )  )
157            DO jk = 2, jpkm1
158               DO jj = 1, jpj
159                  DO ji = 1, jpi
160                     IF( fsdept(ji,jj,jk) <= zh ) THEN
161                        ahm3(ji,jj,jk) = ahm3(ji,jj,1)
162                        ahm4(ji,jj,jk) = ahm4(ji,jj,1)
163                     ELSE
164                        zd = 1.e0 + ( zr - 1.e0 ) * (  1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh )  ) * zc
165                        ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd
166                        ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd
167                     ENDIF
168                  END DO
169               END DO
170            END DO
171            ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
172            ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
173            IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
174            DO jk = 1, jpk
175               IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk
176            END DO
177         ENDIF
178
179         ! Control print
180         IF( lwp .AND. ld_print ) THEN
181            WRITE(numout,*)
182            WRITE(numout,*) 'inildf: ahm3 array at level 1'
183            CALL prihre(ahm3(:,:,1  ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
184            WRITE(numout,*)
185            WRITE(numout,*) 'inildf: ahm4 array at level 1'
186            CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
187         ENDIF
188      ENDIF
189
190      IF(.not. wrk_release(1,2))THEN
191         CALL ctl_stop('ldf_dyn_c3d: ERROR: failed to release workspace array.')
192      END IF
193
194   END SUBROUTINE ldf_dyn_c3d
195
196
197   SUBROUTINE ldf_dyn_c3d_orca( ld_print )
198      !!----------------------------------------------------------------------
199      !!                  ***  ROUTINE ldf_dyn_c3d  ***
200      !!                   
201      !! ** Purpose :   ORCA R1, R2 and R4 only
202      !!
203      !! ** Method  :   blah blah blah ....
204      !!----------------------------------------------------------------------
205      USE ldftra_oce, ONLY : aht0
206      USE wrk_nemo, ONLY: wrk_use, wrk_release, iwrk_use, iwrk_release
207      USE wrk_nemo, ONLY:  icof => iwrk_2d_1
208      USE wrk_nemo, ONLY: zahm0 =>  wrk_2d_1
209      USE wrk_nemo, ONLY: zcoef =>  wrk_1d_1
210      !!
211      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
212      !!
213      INTEGER ::   ji, jj, jk, jn      ! dummy loop indices
214      INTEGER ::   ii0, ii1, ij0, ij1  ! temporary integers
215      INTEGER ::   inum                ! temporary logical unit
216      INTEGER ::   iim, ijm
217      INTEGER ::   ifreq, il1, il2, ij, ii
218      INTEGER, DIMENSION(jpidta, jpjdta) ::   idata
219
220      REAL(wp) ::   &
221         zahmeq, zcoff, zcoft, zmsk,   & ! ???
222         zemax, zemin, zeref, zahmm
223
224      CHARACTER (len=15) ::   clexp
225      !!----------------------------------------------------------------------
226
227      IF( (.not. iwrk_use(2,1)) .OR. (.not. wrk_use(2,1)) .OR. &
228          (.not. wrk_use(1,1)))THEN
229         CALL ctl_stop('ldf_dyn_c3d_orca: ERROR: requested workspace arrays are unavailable.')
230         RETURN
231      END IF
232
233      IF(lwp) WRITE(numout,*)
234      IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient'
235      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
236      IF(lwp) WRITE(numout,*)
237      IF(lwp) WRITE(numout,*) '        orca R1, R2 or R4 ocean model'
238      IF(lwp) WRITE(numout,*) '  reduced in the surface Eq. strip '
239      IF(lwp) WRITE(numout,*)
240
241      ! Read 2d integer array to specify western boundary increase in the
242      ! ===================== equatorial strip (20N-20S) defined at t-points
243
244      CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
245      READ(inum,9101) clexp, iim, ijm
246      READ(inum,'(/)')
247      ifreq = 40
248      il1 = 1
249      DO jn = 1, jpidta/ifreq+1
250         READ(inum,'(/)')
251         il2 = MIN( jpidta, il1+ifreq-1 )
252         READ(inum,9201) ( ii, ji = il1, il2, 5 )
253         READ(inum,'(/)')
254         DO jj = jpjdta, 1, -1
255            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
256         END DO
257         il1 = il1 + ifreq
258      END DO
259     
260      DO jj = 1, nlcj
261         DO ji = 1, nlci
262            icof(ji,jj) = idata( mig(ji), mjg(jj) )
263         END DO
264      END DO
265      DO jj = nlcj+1, jpj
266         DO ji = 1, nlci
267            icof(ji,jj) = icof(ji,nlcj)
268         END DO
269      END DO
270      DO jj = 1, jpj
271         DO ji = nlci+1, jpi
272            icof(ji,jj) = icof(nlci,jj)
273         END DO
274      END DO
275     
2769101  FORMAT(1x,a15,2i8)
2779201  FORMAT(3x,13(i3,12x))
2789202  FORMAT(i3,41i3)
279     
280      ! Set ahm1 and ahm2
281      ! =================
282     
283      ! define ahm1 and ahm2 at the right grid point position
284      ! (USER: modify ahm1 and ahm2 following your desiderata)
285      ! biharmonic : ahm1 (ahm2) defined at u- (v-) point
286      ! harmonic   : ahm1 (ahm2) defined at t- (f-) point
287     
288      ! first level : as for 2D coefficients
289     
290      ! Decrease ahm to zahmeq m2/s in the tropics
291      ! (from 90 to 20 degre: ahm = constant
292      ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
293      ! from  2.5 to  0 degre: ahm = constant
294      ! symmetric in the south hemisphere)
295     
296      IF( jp_cfg == 4 )   THEN
297         zahmeq = 5.0 * aht0
298         zahmm  = min( 160000.0, ahm0)
299         zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
300         zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
301         zeref = MAXVAL ( e1t(:,:) * e2t(:,:),   &
302             &   tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. )
303 
304         DO jj = 1, jpj
305           DO ji = 1, jpi
306              zmsk = e1t(ji,jj) * e2t(ji,jj)
307              IF( abs(gphit(ji,jj)) .LE. 15 ) THEN
308                 zahm0(ji,jj) = ahm0
309              ELSE
310                 IF( zmsk .GE. zeref ) THEN
311                    zahm0(ji,jj) = ahm0
312                 ELSE
313                    zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 -   &
314                        &          cos((rpi*0.5*(zmsk-zemin)/(zeref-zemin))))
315                 ENDIF
316              ENDIF
317           END DO
318         END DO
319      ENDIF
320
321      IF( jp_cfg == 2 )   THEN
322         zahmeq     = aht0
323         zahmm      = ahm0
324         zahm0(:,:) = ahm0
325      ENDIF
326
327      IF( jp_cfg == 1 )   THEN
328         zahmeq     = aht0  ! reduced to aht0 on equator; set to ahm0 if no tropical reduction is required
329         zahmm      = ahm0
330         zahm0(:,:) = ahm0
331      ENDIF
332
333      DO jj = 1, jpj
334         DO ji = 1, jpi
335            IF( ABS(gphif(ji,jj)) >= 20.) THEN
336               ahm2(ji,jj,1) =  zahm0(ji,jj)
337            ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
338               ahm2(ji,jj,1) =  zahmeq
339            ELSE
340               ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2.   &
341                  &            *(1.-COS( rad*(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
342            ENDIF
343            IF( ABS(gphit(ji,jj)) >= 20.) THEN
344               ahm1(ji,jj,1) =  zahm0(ji,jj)
345            ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
346               ahm1(ji,jj,1) =  zahmeq
347            ELSE
348               ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2.   &
349                  &            *(1.-COS( rad*(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
350            ENDIF
351         END DO
352      END DO
353     
354      ! increase along western boundaries of equatorial strip
355      ! t-point
356      DO jj = 1, jpjm1
357         DO ji = 1, jpim1
358            zcoft = float( icof(ji,jj) ) / 100.
359            ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1)
360         END DO
361      END DO
362      ! f-point
363      icof(:,:) = icof(:,:) * tmask(:,:,1)
364      DO jj = 1, jpjm1
365         DO ji = 1, jpim1   ! NO vector opt.
366            zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
367            IF( zmsk == 0. ) THEN
368               zcoff = 1.
369            ELSE
370               zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
371                     / (zmsk * 100.)
372            ENDIF
373            ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1)
374         END DO
375      END DO
376
377      ! other level: re-increase the coef in the deep ocean
378      !==================================================================
379      ! Prior to v3.3, zcoeff was hardwired according to k-index jk.
380      !
381      ! From v3.3 onwards this has been generalised to a function of
382      ! depth so that it can be used with any number of levels.
383      !
384      ! The function has been chosen to match the original values (shown
385      ! in the following comments) when using the standard 31 ORCA levels. 
386      !     DO jk = 1, 21
387      !        zcoef(jk) = 1._wp
388      !     END DO
389      !     zcoef(22) = 2._wp
390      !     zcoef(23) = 3._wp
391      !     zcoef(24) = 5._wp
392      !     zcoef(25) = 7._wp
393      !     zcoef(26) = 9._wp
394      !     DO jk = 27, jpk
395      !        zcoef(jk) = 10._wp
396      !     END DO
397      !==================================================================
398
399       IF(lwp) THEN
400          WRITE(numout,*)
401          WRITE(numout,*) '         1D zcoef array '
402          WRITE(numout,*) '         ~~~~~~~~~~~~~~ '
403          WRITE(numout,*)
404          WRITE(numout,*) '    jk        zcoef '
405       ENDIF
406
407      DO jk=1, jpk
408         zcoef(jk) = 1.0_wp + NINT(9.0_wp*(gdept_0(jk)-800.0_wp)/(3000.0_wp-800.0_wp))
409         zcoef(jk) = MIN(10.0_wp, MAX(1.0_wp, zcoef(jk)))
410         IF(lwp) WRITE(numout,'(4x,i3,6x,f7.3)') jk,zcoef(jk)
411      END DO
412
413      DO jk = 2, jpk
414         ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) )
415         ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) )
416      END DO
417
418      IF( jp_cfg == 4 )   THEN         ! Limit AHM in Gibraltar strait
419         ij0 = 50   ;   ij1 = 53
420         ii0 = 69   ;   ii1 = 71
421         DO jk = 1, jpk
422            ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
423            ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
424         END DO
425      ENDIF
426      CALL lbc_lnk( ahm1, 'T', 1. )   ! Lateral boundary conditions (unchanged sign)
427      CALL lbc_lnk( ahm2, 'F', 1. )
428
429
430      IF(lwp) THEN                    ! Control print
431         WRITE(numout,*)
432         WRITE(numout,*) '         3D ahm1 array (k=1)'
433         CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
434         WRITE(numout,*)
435         WRITE(numout,*) '         3D ahm2 array (k=1)'
436         CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
437         WRITE(numout,*)
438         WRITE(numout,*) '         3D ahm2 array (k=jpk)'
439         CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
440      ENDIF
441
442
443      ! Set ahm3 and ahm4
444      ! =================
445
446      ! define ahm3 and ahm4 at the right grid point position
447      ! initialization to a constant value
448      !     (USER: modify ahm3 and ahm4 following your desiderata)
449      !     harmonic isopycnal or geopotential:
450      !                          ahm3 (ahm4) defined at u- (v-) point
451      DO jk = 1, jpk
452         DO jj = 2, jpj
453            DO ji = 2, jpi
454               ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji  ,jj-1,jk) )
455               ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj  ,jk) )
456            END DO
457         END DO
458      END DO
459      ahm3 ( :, 1, :) = ahm3 ( :, 2, :)
460      ahm4 ( :, 1, :) = ahm4 ( :, 2, :)
461     
462      CALL lbc_lnk( ahm3, 'U', 1. )    ! Lateral boundary conditions (unchanged sign)
463      CALL lbc_lnk( ahm4, 'V', 1. )
464
465      ! Control print
466
467      IF( lwp .AND. ld_print ) THEN
468         WRITE(numout,*)
469         WRITE(numout,*) '         ahm3 array level 1'
470         CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
471         WRITE(numout,*)
472         WRITE(numout,*) '         ahm4 array level 1'
473         CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
474      ENDIF
475
476      IF( (.not. iwrk_release(2,1)) .OR. (.not. wrk_release(2,1)) .OR. &
477          (.not. wrk_release(1,1)))THEN
478         CALL ctl_stop('ldf_dyn_c3d_orca: ERROR: failed to release workspace arrays.')
479      END IF
480
481   END SUBROUTINE ldf_dyn_c3d_orca
Note: See TracBrowser for help on using the repository browser.