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.
domvvl.F90 in branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 2828

Last change on this file since 2828 was 2828, checked in by acc, 13 years ago

Branch: dev_r2802_NOCS_vvlfix. Bugfix to vvl code when using modified metrics in key straits (ORCA 2, 1 and 05 configurations). Calculations of e3u,v metrics in domvvl_2 now ignore any modified metrics when area averaging. This fixes castrophic errors that were occurring at the modified Straits. The energetic consistency is already broken by the modification of the metric and this change does not introduce any addition error.
  • Property svn:keywords set to Id
File size: 22.4 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!----------------------------------------------------------------------
9#if defined key_vvl
10   !!----------------------------------------------------------------------
11   !!   'key_vvl'                              variable volume
12   !!----------------------------------------------------------------------
13   !!   dom_vvl     : defined coefficients to distribute ssh on each layers
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE sbc_oce         ! surface boundary condition: ocean
18   USE phycst          ! physical constants
19   USE in_out_manager  ! I/O manager
20   USE lib_mpp         ! distributed memory computing library
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   dom_vvl         ! called by domain.F90
27   PUBLIC   dom_vvl_2       ! called by domain.F90
28   PUBLIC   dom_vvl_alloc   ! called by nemogcm.F90
29
30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: 1/H_0 at t-,u-,v-,f-points
31
32   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra
33      !                                                              ! except at nit000 (=rdttra) if neuler=0
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS       
44
45   INTEGER FUNCTION dom_vvl_alloc()
46      !!----------------------------------------------------------------------
47      !!                ***  ROUTINE dom_vvl_alloc  ***
48      !!----------------------------------------------------------------------
49      !
50      ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) ,     &
51         &      r2dt        (jpk)                                                             , STAT=dom_vvl_alloc )
52         !
53      IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
54      IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
55      !
56   END FUNCTION dom_vvl_alloc
57
58
59   SUBROUTINE dom_vvl
60      !!----------------------------------------------------------------------
61      !!                ***  ROUTINE dom_vvl  ***
62      !!                   
63      !! ** Purpose :   compute mu coefficients at t-, u-, v- and f-points to
64      !!              spread ssh over the whole water column (scale factors)
65      !!                set the before and now ssh at u- and v-points
66      !!              (also f-point in now case)
67      !!----------------------------------------------------------------------
68      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
69      USE wrk_nemo, ONLY:   zee_t => wrk_2d_1, zee_u => wrk_2d_2, zee_v => wrk_2d_3, zee_f => wrk_2d_4   ! 2D workspace
70      !
71      INTEGER  ::   ji, jj, jk   ! dummy loop indices
72      REAL(wp) ::   zcoefu, zcoefv , zcoeff                ! local scalars
73      REAL(wp) ::   zvt   , zvt_ip1, zvt_jp1, zvt_ip1jp1   !   -      -
74      !!----------------------------------------------------------------------
75
76      IF( wrk_in_use(2, 1,2,3,4) ) THEN
77         CALL ctl_stop('dom_vvl: requested workspace arrays unavailable')   ;   RETURN
78      ENDIF
79
80      IF(lwp) THEN
81         WRITE(numout,*)
82         WRITE(numout,*) 'dom_vvl : Variable volume initialization'
83         WRITE(numout,*) '~~~~~~~~  compute coef. used to spread ssh over each layers'
84      ENDIF
85     
86      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl : unable to allocate arrays' )
87
88      fsdept(:,:,:) = gdept (:,:,:)
89      fsdepw(:,:,:) = gdepw (:,:,:)
90      fsde3w(:,:,:) = gdep3w(:,:,:)
91      fse3t (:,:,:) = e3t   (:,:,:)
92      fse3u (:,:,:) = e3u   (:,:,:)
93      fse3v (:,:,:) = e3v   (:,:,:)
94      fse3f (:,:,:) = e3f   (:,:,:)
95      fse3w (:,:,:) = e3w   (:,:,:)
96      fse3uw(:,:,:) = e3uw  (:,:,:)
97      fse3vw(:,:,:) = e3vw  (:,:,:)
98
99      !                                 !==  mu computation  ==!
100      zee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level
101      zee_u(:,:) = fse3u_0(:,:,1)
102      zee_v(:,:) = fse3v_0(:,:,1)
103      zee_f(:,:) = fse3f_0(:,:,1)
104      DO jk = 2, jpkm1                          ! Sum of the masked vertical scale factors
105         zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
106         zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
107         zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
108         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask
109            zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk)
110         END DO
111      END DO 
112      !                                         ! Compute and mask the inverse of the local depth at T, U, V and F points
113      zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1)
114      zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1)
115      zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1)
116      DO jj = 1, jpjm1                               ! f-point case fmask cannot be used
117         zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1)
118      END DO
119      CALL lbc_lnk( zee_f, 'F', 1. )                 ! lateral boundary condition on ee_f
120      !
121      DO jk = 1, jpk                            ! mu coefficients
122         mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels
123         muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk)     ! U-point at T levels
124         muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels
125      END DO
126      DO jk = 1, jpk                                 ! F-point : fmask=shlat at coasts, use the product of umask
127         DO jj = 1, jpjm1
128               muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels
129         END DO
130         muf(:,jpj,jk) = 0._wp
131      END DO
132      CALL lbc_lnk( muf, 'F', 1. )                   ! lateral boundary condition
133
134
135      hu_0(:,:) = 0.e0                          ! Reference ocean depth at U- and V-points
136      hv_0(:,:) = 0.e0
137      DO jk = 1, jpk
138         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
139         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
140      END DO
141     
142      DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
143         DO ji = 1, jpim1   ! NO vector opt.
144            zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1)
145            zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1)
146            zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1)
147            !
148            zvt           = e1e2t(ji  ,jj  ) * sshb(ji  ,jj  )    ! before fields
149            zvt_ip1       = e1e2t(ji+1,jj  ) * sshb(ji+1,jj  )
150            zvt_jp1       = e1e2t(ji  ,jj+1) * sshb(ji  ,jj+1)
151            sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 )
152            sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 )
153            !
154            zvt           = e1e2t(ji  ,jj  ) * sshn(ji  ,jj  )    ! now fields
155            zvt_ip1       = e1e2t(ji+1,jj  ) * sshn(ji+1,jj  )
156            zvt_jp1       = e1e2t(ji  ,jj+1) * sshn(ji  ,jj+1)
157            zvt_ip1jp1    = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1)
158            sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 )
159            sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 )
160            sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 )
161         END DO
162      END DO
163      CALL lbc_lnk( sshu_n, 'U', 1. )   ;   CALL lbc_lnk( sshu_b, 'U', 1. )      ! lateral boundary conditions
164      CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. )
165      CALL lbc_lnk( sshf_n, 'F', 1. )
166      !
167      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays')
168      !
169   END SUBROUTINE dom_vvl
170
171
172   SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b )
173      !!----------------------------------------------------------------------
174      !!                ***  ROUTINE dom_vvl_2  ***
175      !!                   
176      !! ** Purpose :   compute the vertical scale factors at u- and v-points
177      !!              in variable volume case.
178      !!
179      !! ** Method  :   In variable volume case (non linear sea surface) the
180      !!              the vertical scale factor at velocity points is computed
181      !!              as the average of the cell surface weighted e3t.
182      !!                It uses the sea surface heigth so it have to be initialized
183      !!              after ssh is read/set
184      !!----------------------------------------------------------------------
185      INTEGER                   , INTENT(in   ) ::   kt               ! ocean time-step index
186      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3u_b, pe3v_b   ! before vertical scale factor at u- & v-pts
187      !
188      INTEGER  ::   ji, jj, jk   ! dummy loop indices
189      INTEGER  ::   iku, ikv     ! local integers   
190      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
191      REAL(wp) ::   zvt          ! local scalars
192      !!----------------------------------------------------------------------
193
194      IF( lwp .AND. kt == nit000 ) THEN
195         WRITE(numout,*)
196         WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization'
197         WRITE(numout,*) '~~~~~~~~~ '
198         pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk)
199         pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk)
200      ENDIF
201     
202      DO jk = 1, jpkm1           ! set the before scale factors at u- & v-points
203         DO jj = 2, jpjm1
204            DO ji = fs_2, fs_jpim1
205               zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj)
206               pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) )
207               pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) )
208            END DO
209         END DO
210      END DO
211
212      ! Correct scale factors at locations that have been individually modified in domhgr
213      ! Such modifications break the relationship between e1e2t and e1u*e2u etc. Recompute
214      ! scale factors ignoring the modified metric.
215      !                                                ! =====================
216      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
217         !                                             ! =====================
218         IF( nn_cla == 0 ) THEN
219            !
220            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
221            ij0 = 102   ;   ij1 = 102   
222            DO jk = 1, jpkm1                 ! set the before scale factors at u-points
223               DO jj = mj0(ij0), mj1(ij1)
224                  DO ji = mi0(ii0), mi1(ii1)
225                     zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
226                     pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
227                  END DO
228               END DO
229            END DO
230            !
231            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
232            ij0 =  88   ;   ij1 =  88   
233            DO jk = 1, jpkm1                 ! set the before scale factors at u-points
234               DO jj = mj0(ij0), mj1(ij1)
235                  DO ji = mi0(ii0), mi1(ii1)
236                     zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
237                     pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
238                  END DO
239               END DO
240            END DO
241            DO jk = 1, jpkm1                 ! set the before scale factors at v-points
242               DO jj = mj0(ij0), mj1(ij1)
243                  DO ji = mi0(ii0), mi1(ii1)
244                     zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
245                     pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
246                  END DO
247               END DO
248            END DO
249         ENDIF
250
251         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
252         ij0 = 116   ;   ij1 = 116   
253         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
254            DO jj = mj0(ij0), mj1(ij1)
255               DO ji = mi0(ii0), mi1(ii1)
256                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
257                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
258               END DO
259            END DO
260         END DO
261         !
262      ENDIF
263         !                                             ! =====================
264      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
265         !                                             ! =====================
266
267         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
268         ij0 = 200   ;   ij1 = 200   
269         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
270            DO jj = mj0(ij0), mj1(ij1)
271               DO ji = mi0(ii0), mi1(ii1)
272                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
273                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
274               END DO
275            END DO
276         END DO
277
278         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
279         ij0 = 208   ;   ij1 = 208   
280         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
281            DO jj = mj0(ij0), mj1(ij1)
282               DO ji = mi0(ii0), mi1(ii1)
283                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
284                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
285               END DO
286            END DO
287         END DO
288
289         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
290         ij0 = 124   ;   ij1 = 125   
291         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
292            DO jj = mj0(ij0), mj1(ij1)
293               DO ji = mi0(ii0), mi1(ii1)
294                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
295                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
296               END DO
297            END DO
298         END DO
299
300         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
301         ij0 = 124   ;   ij1 = 125   
302         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
303            DO jj = mj0(ij0), mj1(ij1)
304               DO ji = mi0(ii0), mi1(ii1)
305                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
306                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
307               END DO
308            END DO
309         END DO
310
311         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
312         ij0 = 124   ;   ij1 = 125   
313         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
314            DO jj = mj0(ij0), mj1(ij1)
315               DO ji = mi0(ii0), mi1(ii1)
316                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
317                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
318               END DO
319            END DO
320         END DO
321
322         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
323         ij0 = 124   ;   ij1 = 125   
324         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
325            DO jj = mj0(ij0), mj1(ij1)
326               DO ji = mi0(ii0), mi1(ii1)
327                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
328                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
329               END DO
330            END DO
331         END DO
332
333         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
334         ij0 = 141   ;   ij1 = 142   
335         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
336            DO jj = mj0(ij0), mj1(ij1)
337               DO ji = mi0(ii0), mi1(ii1)
338                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
339                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
340               END DO
341            END DO
342         END DO
343
344         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
345         ij0 = 141   ;   ij1 = 142   
346         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
347            DO jj = mj0(ij0), mj1(ij1)
348               DO ji = mi0(ii0), mi1(ii1)
349                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
350                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
351               END DO
352            END DO
353         END DO
354
355         !
356      ENDIF
357      !                                                ! ======================
358      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
359         !                                             ! ======================
360         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
361         ij0 = 327   ;   ij1 = 327   
362         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
363            DO jj = mj0(ij0), mj1(ij1)
364               DO ji = mi0(ii0), mi1(ii1)
365                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
366                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
367               END DO
368            END DO
369         END DO
370         !
371         ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u was modified)
372         ij0 = 343   ;   ij1 = 343   
373         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
374            DO jj = mj0(ij0), mj1(ij1)
375               DO ji = mi0(ii0), mi1(ii1)
376                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
377                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
378               END DO
379            END DO
380         END DO
381         !
382         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
383         ij0 = 232   ;   ij1 = 232   
384         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
385            DO jj = mj0(ij0), mj1(ij1)
386               DO ji = mi0(ii0), mi1(ii1)
387                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
388                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
389               END DO
390            END DO
391         END DO
392         !
393         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
394         ij0 = 232   ;   ij1 = 232   
395         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
396            DO jj = mj0(ij0), mj1(ij1)
397               DO ji = mi0(ii0), mi1(ii1)
398                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
399                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
400               END DO
401            END DO
402         END DO
403         !
404         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
405         ij0 = 270   ;   ij1 = 270   
406         DO jk = 1, jpkm1                 ! set the before scale factors at u-points
407            DO jj = mj0(ij0), mj1(ij1)
408               DO ji = mi0(ii0), mi1(ii1)
409                  zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj)
410                  pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) )
411               END DO
412            END DO
413         END DO
414         !
415         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
416         ij0 = 232   ;   ij1 = 233   
417         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
418            DO jj = mj0(ij0), mj1(ij1)
419               DO ji = mi0(ii0), mi1(ii1)
420                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
421                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
422               END DO
423            END DO
424         END DO
425         !
426         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
427         ij0 = 276   ;   ij1 = 276   
428         DO jk = 1, jpkm1                 ! set the before scale factors at v-points
429            DO jj = mj0(ij0), mj1(ij1)
430               DO ji = mi0(ii0), mi1(ii1)
431                  zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)
432                  pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )
433               END DO
434            END DO
435         END DO
436         !
437      ENDIF
438      ! End of individual corrections to scale factors
439
440      IF( ln_zps ) THEN          ! minimum of the e3t at partial cell level
441         DO jj = 2, jpjm1
442            DO ji = fs_2, fs_jpim1
443               iku = mbku(ji,jj)
444               ikv = mbkv(ji,jj)
445               pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj  ,iku) ) 
446               pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji  ,jj+1,ikv) ) 
447            END DO
448         END DO
449      ENDIF
450
451      pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:)      ! anomaly to avoid zero along closed boundary/extra halos
452      pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:)
453      CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
454      CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. )
455      pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:)      ! recover the full scale factor
456      pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:)
457      !
458   END SUBROUTINE dom_vvl_2
459   
460#else
461   !!----------------------------------------------------------------------
462   !!   Default option :                                      Empty routine
463   !!----------------------------------------------------------------------
464CONTAINS
465   SUBROUTINE dom_vvl
466   END SUBROUTINE dom_vvl
467   SUBROUTINE dom_vvl_2(kdum, pudum, pvdum )
468      USE par_kind
469      INTEGER                   , INTENT(in   ) ::   kdum       
470      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pudum, pvdum
471   END SUBROUTINE dom_vvl_2
472#endif
473
474   !!======================================================================
475END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.