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.
cla.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LBC – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LBC/cla.F90 @ 2636

Last change on this file since 2636 was 2636, checked in by gm, 13 years ago

dynamic mem: #785 ; move ctl_stop & warn in lib_mpp to avoid a circular dependency + ctl_stop improvment

  • Property svn:keywords set to Id
File size: 41.6 KB
Line 
1MODULE cla
2   !!======================================================================
3   !!                    ***  MODULE  cla  ***
4   !! Cross Land Advection : specific update of the horizontal divergence,
5   !!                        tracer trends and after velocity
6   !!
7   !!                 ---   Specific to ORCA_R2   ---
8   !!
9   !!======================================================================
10   !! History :  1.0  ! 2002-11 (A. Bozec)  Original code
11   !!            3.2  ! 2009-07 (G. Madec)  merge cla, cla_div, tra_cla, cla_dynspg
12   !!                 !                     and correct a mpp bug reported by A.R. Porter
13   !!----------------------------------------------------------------------
14#if defined key_orca_r2
15   !!----------------------------------------------------------------------
16   !!   'key_orca_r2'                                 global ocean model R2
17   !!----------------------------------------------------------------------
18   !!   cla_div           : update of horizontal divergence at cla straits
19   !!   tra_cla           : update of tracers at cla straits
20   !!   cla_dynspg        : update of after horizontal velocities at cla straits
21   !!   cla_init          : initialisation - control check
22   !!   cla_bab_el_mandeb : cross land advection for Bab-el-mandeb strait
23   !!   cla_gibraltar     : cross land advection for Gibraltar strait
24   !!   cla_hormuz        : cross land advection for Hormuz strait
25   !!----------------------------------------------------------------------
26   USE oce            ! ocean dynamics and tracers
27   USE dom_oce        ! ocean space and time domain
28   USE sbc_oce        ! surface boundary condition: ocean
29   USE dynspg_oce     ! ocean dynamics: surface pressure gradient variables
30   USE in_out_manager ! I/O manager
31   USE lib_mpp        ! distributed memory computing library
32   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
33   USE lib_mpp        ! MPP library
34
35   IMPLICIT NONE
36   PRIVATE
37   
38   PUBLIC   cla_init     ! routine called by opa.F90
39   PUBLIC   cla_div      ! routine called by divcur.F90
40   PUBLIC   cla_traadv   ! routine called by traadv.F90
41   PUBLIC   cla_dynspg   ! routine called by dynspg_flt.F90
42
43   INTEGER  ::   nbab, ngib, nhor   ! presence or not of required grid-point on local domain
44   !                                ! for Bab-el-Mandeb, Gibraltar, and Hormuz straits
45   
46   !                                                              !!! profile of hdiv for some straits
47   REAL(wp), DIMENSION (jpk) ::   hdiv_139_101, hdiv_139_101_kt    ! Gibraltar     strait, fixed & time evolving part (i,j)=(172,101)
48   REAL(wp), DIMENSION (jpk) ::   hdiv_139_102                     ! Gibraltar     strait, fixed part only            (i,j)=(139,102)
49   REAL(wp), DIMENSION (jpk) ::   hdiv_141_102, hdiv_141_102_kt    ! Gibraltar     strait, fixed & time evolving part (i,j)=(141,102)
50   REAL(wp), DIMENSION (jpk) ::   hdiv_161_88 , hdiv_161_88_kt     ! Bab-el-Mandeb strait, fixed & time evolving part (i,j)=(161,88)
51   REAL(wp), DIMENSION (jpk) ::   hdiv_161_87                      ! Bab-el-Mandeb strait, fixed part only            (i,j)=(161,87)
52   REAL(wp), DIMENSION (jpk) ::   hdiv_160_89 , hdiv_160_89_kt     ! Bab-el-Mandeb strait, fixed & time evolving part (i,j)=(160,89)
53   REAL(wp), DIMENSION (jpk) ::   hdiv_172_94                      ! Hormuz        strait, fixed part only            (i,j)=(172, 94)
54
55   REAL(wp), DIMENSION (jpk) ::   t_171_94_hor, s_171_94_hor       ! Temperature, salinity in the Hormuz strait
56   
57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE cla_div( kt )
67      !!----------------------------------------------------------------------
68      !!                 ***  ROUTINE div_cla  ***
69      !!
70      !! ** Purpose :   update the horizontal divergence of the velocity field
71      !!              at some straits ( Gibraltar, Bab el Mandeb and Hormuz ).
72      !!
73      !! ** Method  : - first time-step: initialisation of cla
74      !!              - all   time-step: using imposed transport at each strait,
75      !!              the now horizontal divergence is updated
76      !!
77      !! ** Action  :   phdivn   updted now horizontal divergence at cla straits
78      !!----------------------------------------------------------------------
79      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
80      !!----------------------------------------------------------------------
81      !     
82      IF( kt == nit000 ) THEN
83         !
84         CALL cla_init                                        ! control check
85         !
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'div_cla : cross land advection on hdiv '
88         IF(lwp) WRITE(numout,*) '~~~~~~~~'
89         !
90         IF( nbab == 1 )   CALL cla_bab_el_mandeb('ini')    ! Bab el Mandeb ( Red Sea - Indian ocean )
91         IF( ngib == 1 )   CALL cla_gibraltar    ('ini')    ! Gibraltar strait (Med Sea - Atlantic ocean)
92         IF( nhor == 1 )   CALL cla_hormuz       ('ini')    ! Hormuz Strait ( Persian Gulf - Indian ocean )
93         !
94      ENDIF                           
95      !
96      IF( nbab == 1    )   CALL cla_bab_el_mandeb('div')    ! Bab el Mandeb ( Red Sea - Indian ocean )
97      IF( ngib == 1    )   CALL cla_gibraltar    ('div')    ! Gibraltar strait (Med Sea - Atlantic ocean)
98      IF( nhor == 1    )   CALL cla_hormuz       ('div')    ! Hormuz Strait ( Persian Gulf - Indian ocean )
99      !
100!!gm  lbc useless here, no?
101!!gm      CALL lbc_lnk( hdivn, 'T', 1. )                    ! Lateral boundary conditions on hdivn
102      !
103   END SUBROUTINE cla_div
104   
105   
106   SUBROUTINE cla_traadv( kt )
107      !!----------------------------------------------------------------------
108      !!                 ***  ROUTINE tra_cla  ***
109      !!                   
110      !! ** Purpose :   Update the now trend due to the advection of tracers
111      !!      and add it to the general trend of passive tracer equations
112      !!      at some straits ( Bab el Mandeb, Gibraltar, Hormuz ).
113      !!
114      !! ** Method  :   using both imposed transport at each strait and T & S
115      !!              budget, the now tracer trends is updated
116      !!
117      !! ** Action  :   (ta,sa)   updated now tracer trends at cla straits
118      !!----------------------------------------------------------------------
119      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
120      !!----------------------------------------------------------------------
121      !
122      IF( kt == nit000 ) THEN
123         IF(lwp) WRITE(numout,*)
124         IF(lwp) WRITE(numout,*) 'tra_cla : cross land advection on tracers '
125         IF(lwp) WRITE(numout,*) '~~~~~~~~'
126      ENDIF
127      !
128      IF( nbab == 1    )   CALL cla_bab_el_mandeb('tra')      ! Bab el Mandeb strait
129      IF( ngib == 1    )   CALL cla_gibraltar    ('tra')      ! Gibraltar strait
130      IF( nhor == 1    )   CALL cla_hormuz       ('tra')      ! Hormuz Strait ( Persian Gulf)
131      !
132   END SUBROUTINE cla_traadv
133
134   
135   SUBROUTINE cla_dynspg( kt )
136      !!----------------------------------------------------------------------
137      !!                 ***  ROUTINE cla_dynspg  ***
138      !!                   
139      !! ** Purpose :   Update the after velocity at some straits
140      !!              (Bab el Mandeb, Gibraltar, Hormuz).
141      !!
142      !! ** Method  :   required to compute the filtered surface pressure gradient
143      !!
144      !! ** Action  :   (ua,va)   after velocity at the cla straits
145      !!----------------------------------------------------------------------
146      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
147      !!----------------------------------------------------------------------
148      !
149      IF( kt == nit000 ) THEN
150         IF(lwp) WRITE(numout,*)
151         IF(lwp) WRITE(numout,*) 'cla_dynspg : cross land advection on (ua,va) '
152         IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
153      ENDIF
154      !
155      IF( nbab == 1    )   CALL cla_bab_el_mandeb('spg')      ! Bab el Mandeb strait
156      IF( ngib == 1    )   CALL cla_gibraltar    ('spg')      ! Gibraltar strait
157      IF( nhor == 1    )   CALL cla_hormuz       ('spg')      ! Hormuz Strait ( Persian Gulf)
158      !
159!!gm lbc is needed here, not?
160!!gm      CALL lbc_lnk( hdivn, 'U', -1. )   ;   CALL lbc_lnk( hdivn, 'V', -1. )      ! Lateral boundary conditions
161      !
162   END SUBROUTINE cla_dynspg
163
164
165   SUBROUTINE cla_init
166      !! -------------------------------------------------------------------
167      !!                   ***  ROUTINE cla_init  ***
168      !!           
169      !! ** Purpose :   control check for mpp computation 
170      !!
171      !! ** Method  : - All the strait grid-points must be inside one of the
172      !!              local domain interior for the cla advection to work
173      !!              properly in mpp (i.e. inside (2:jpim1,2:jpjm1) ).
174      !!              Define the corresponding indicators (nbab, ngib, nhor)
175      !!              - The profiles of cross-land fluxes are currently hard
176      !!              coded for L31 levels. Stop if jpk/=31
177      !!
178      !! ** Action  :   nbab, ngib, nhor   strait inside the local domain or not
179      !!---------------------------------------------------------------------
180      REAL(wp) ::   ztemp
181      !!---------------------------------------------------------------------
182      !
183      IF(lwp) WRITE(numout,*)
184      IF(lwp) WRITE(numout,*) 'cla_init : cross land advection initialisation '
185      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
186      !
187      IF( .NOT.lk_dynspg_flt )   CALL ctl_stop( 'cla_init: Cross Land Advection works only with lk_dynspg_flt=T ' )
188      !
189      IF( lk_vvl    )   CALL ctl_stop( 'cla_init: Cross Land Advection does not work with lk_vvl=T option' )
190      !
191      IF( jpk /= 31 )   CALL ctl_stop( 'cla_init: Cross Land Advection hard coded for ORCA_R2_L31' )
192      !
193      !                                        _|_______|_______|_
194      !                                     89  |       |///////| 
195      !                                        _|_______|_______|_
196      ! ------------------------ !          88  |///////|       |
197      !   Bab el Mandeb strait   !             _|_______|_______|_
198      ! ------------------------ !          87  |///////|       |
199      !                                        _|_______|_______|_
200      !                                         |  160  |  161  | 
201      !
202      ! The 6 Bab el Mandeb grid-points must be inside one of the interior of the
203      ! local domain for the cla advection to work properly (i.e. (2:jpim1,2:jpjm1)
204      nbab = 0
205      IF(  ( 1 <= mj0( 88) .AND. mj1( 89) <= jpj ) .AND.    &  !* (161,89), (161,88) and (161,88) on the local pocessor
206         & ( 1 <= mi0(160) .AND. mi1(161) <= jpi )       )    nbab = 1 
207      !
208      ! test if there is no local domain that includes all required grid-points
209      ztemp = REAL( nbab )
210      IF( lk_mpp )   CALL mpp_sum( ztemp )      ! sum with other processors value
211      IF( ztemp == 0 ) THEN                     ! Only 2 points in each direction, this should never be a problem
212         CALL ctl_stop( ' cross land advection at Bab-el_Mandeb does not work with your processor cutting: change it' )
213      ENDIF
214      !                                        ___________________________
215      ! ------------------------ !         102  |       |///////|       |
216      !     Gibraltar strait     !             _|_______|_______|_______|_
217      ! ------------------------ !         101  |       |///////|       |
218      !                                        _|_______|_______|_______|_
219      !                                         |  139  |  140  |  141  |
220      !
221      ! The 6 Gibraltar grid-points must be inside one of the interior of the
222      ! local domain for the cla advection to work properly (i.e. (2:jpim1,2:jpjm1)
223      ngib = 0
224      IF(  ( 2 <= mj0(101) .AND. mj1(102) <= jpjm1 ) .AND.    &  !* (139:141,101:102) on the local pocessor
225         & ( 2 <= mi0(139) .AND. mi1(141) <= jpim1 )       )    ngib = 1 
226      !
227      ! test if there is no local domain that includes all required grid-points
228      ztemp = REAL( ngib )
229      IF( lk_mpp )   CALL mpp_sum( ztemp )      ! sum with other processors value
230      IF( ztemp == 0 ) THEN                     ! 3 points in i-direction, this may be a problem with some cutting
231           CALL ctl_stop( ' cross land advection at Gibraltar does not work with your processor cutting: change it' )
232      ENDIF
233      !                                        _______________
234      ! ------------------------ !          94  |/////|     |
235      !       Hormuz strait      !             _|_____|_____|_
236      ! ------------------------ !                171   172     
237      !           
238      ! The 2 Hormuz grid-points must be inside one of the interior of the
239      ! local domain for the cla advection to work properly (i.e. (2:jpim1,2:jpjm1)
240      nhor = 0
241      IF(    2 <= mj0( 94) .AND. mj1( 94) <= jpjm1  .AND.  & 
242         &   2 <= mi0(171) .AND. mi1(172) <= jpim1         )   nhor = 1 
243      !
244      ! test if there is no local domain that includes all required grid-points
245      ztemp = REAL( nhor )
246      IF( lk_mpp )   CALL mpp_sum( ztemp )      ! sum with other processors value
247      IF( ztemp == 0 ) THEN                     ! 3 points in i-direction, this may be a problem with some cutting
248           CALL ctl_stop( ' cross land advection at Hormuz does not work with your processor cutting: change it' )
249      ENDIF
250      !
251   END SUBROUTINE cla_init
252
253
254   SUBROUTINE cla_bab_el_mandeb( cd_td )
255      !!----------------------------------------------------------------------
256      !!                ***  ROUTINE cla_bab_el_mandeb  ***
257      !!       
258      !! ** Purpose :   update the now horizontal divergence, the tracer tendancy
259      !!              and the after velocity in vicinity of Bab el Mandeb ( Red Sea - Indian ocean).
260      !!
261      !! ** Method  :   compute the exchanges at each side of the strait :
262      !!
263      !!       surf. zio_flow
264      !! (+ balance of emp) /\  |\\\\\\\\\\\|
265      !!                    ||  |\\\\\\\\\\\| 
266      !!    deep zio_flow   ||  |\\\\\\\\\\\| 
267      !!            |  ||   ||  |\\\\\\\\\\\| 
268      !!        89  |  ||   ||  |\\\\\\\\\\\| 
269      !!            |__\/_v_||__|____________
270      !!            !\\\\\\\\\\\|          surf. zio_flow
271      !!            |\\\\\\\\\\\|<===    (+ balance of emp)
272      !!            |\\\\\\\\\\\u
273      !!        88  |\\\\\\\\\\\|<---      deep  zrecirc (upper+deep at 2 different levels)
274      !!            |___________|__________   
275      !!            !\\\\\\\\\\\|         
276      !!            |\\\\\\\\\\\| ---\     deep  zrecirc (upper+deep)
277      !!        87  !\\\\\\\\\\\u ===/   + deep  zio_flow   (all at the same level)
278      !!            !\\\\\\\\\\\| 
279      !!            !___________|__________
280      !!                160         161
281      !!
282      !!----------------------------------------------------------------------
283      CHARACTER(len=1), INTENT(in) ::   cd_td   ! ='div' update the divergence
284      !                                         ! ='tra' update the tracers
285      !                                         ! ='spg' update after velocity
286      INTEGER  ::   ji, jj, jk   ! dummy loop indices
287      REAL(wp) ::   zemp_red     ! temporary scalar
288      REAL(wp) ::   zio_flow, zrecirc_upp, zrecirc_mid, zrecirc_bot
289      !!---------------------------------------------------------------------
290      !
291      SELECT CASE( cd_td ) 
292      !                     ! ---------------- !
293      CASE( 'ini' )         !  initialisation  !
294         !                  ! ---------------- !
295         !                                   
296         zio_flow    = 0.4e6                       ! imposed in/out flow
297         zrecirc_upp = 0.2e6                       ! imposed upper recirculation water
298         zrecirc_bot = 0.5e6                       ! imposed bottom  recirculation water
299
300         hdiv_161_88(:) = 0.e0                     ! (161,88) Gulf of Aden side, north point
301         hdiv_161_87(:) = 0.e0                     ! (161,87) Gulf of Aden side, south point
302         hdiv_160_89(:) = 0.e0                     ! (160,89) Red sea side
303
304         DO jj = mj0(88), mj1(88)              !** profile of hdiv at (161,88)   (Gulf of Aden side, north point)
305            DO ji = mi0(161), mi1(161)         !------------------------------
306               DO jk = 1, 8                        ! surface in/out flow   (Ind -> Red)   (div >0)
307                  hdiv_161_88(jk) = + zio_flow / ( 8. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
308               END DO
309               !                                   ! recirculation water   (Ind -> Red)   (div >0)
310               hdiv_161_88(20) =                 + zrecirc_upp   / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,20) )
311               hdiv_161_88(21) = + ( zrecirc_bot - zrecirc_upp ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,21) )
312            END DO
313         END DO
314         !
315         DO jj = mj0(87), mj1(87)              !** profile of hdiv at (161,88)   (Gulf of Aden side, south point)
316            DO ji = mi0(161), mi1(161)         !------------------------------
317               !                                   ! deep out flow + recirculation   (Red -> Ind)   (div <0)
318               hdiv_161_87(21) = - ( zio_flow + zrecirc_bot ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,21) )
319            END DO
320         END DO
321         !
322         DO jj = mj0(89), mj1(89)              !** profile of hdiv at (161,88)   (Red sea side)
323            DO ji = mi0(160), mi1(160)         !------------------------------
324               DO jk = 1, 8                        ! surface inflow    (Ind -> Red)   (div <0)
325                  hdiv_160_89(jk) = - zio_flow / ( 8. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
326               END DO
327               !                                   ! deep    outflow   (Red -> Ind)   (div >0)
328               hdiv_160_89(16)    = + zio_flow / (      e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,16) )
329            END DO
330         END DO
331         !                  ! ---------------- !
332      CASE( 'div' )         !   update hdivn   ! (call by divcur module)
333         !                  ! ---------=====-- !
334         !                                     !** emp on the Red Sea   (div >0)
335         zemp_red = 0.e0                       !---------------------
336         DO jj = mj0(87), mj1(96)                  ! sum over the Red sea
337            DO ji = mi0(148), mi1(160) 
338               zemp_red = zemp_red + emp(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
339            END DO
340         END DO
341         IF( lk_mpp )   CALL mpp_sum( zemp_red )   ! sum with other processors value
342         zemp_red = zemp_red * 1.e-3               ! convert in m3
343         !
344         !                                     !** Correct hdivn (including emp adjustment)
345         !                                     !-------------------------------------------
346         DO jj = mj0(88), mj1(88)                  !* profile of hdiv at (161,88)   (Gulf of Aden side, north point)
347            DO ji = mi0(161), mi1(161) 
348               hdiv_161_88_kt(:) = hdiv_161_88(:)
349               DO jk = 1, 8                              ! increase the inflow from the Indian   (div >0)
350                  hdiv_161_88_kt(jk) = hdiv_161_88(jk) + zemp_red / (8. * e2u(ji,jj) * fse3u(ji,jj,jk) )
351               END DO
352               hdivn(ji,jj,:) = hdivn(ji,jj,:) + hdiv_161_88_kt(:)
353            END DO
354         END DO
355         DO jj = mj0(87), mj1(87)                  !* profile of divergence at (161,87)   (Gulf of Aden side, south point)
356            DO ji = mi0(161), mi1(161) 
357               hdivn(ji,jj,:) = hdivn(ji,jj,:) + hdiv_161_87(:)
358            END DO
359         END DO
360         DO jj = mj0(89), mj1(89)                  !* profile of divergence at (160,89)   (Red sea side)
361            DO ji = mi0(160), mi1(160) 
362               hdiv_160_89_kt(:) = hdiv_160_89(:)
363               DO jk = 1, 18                              ! increase the inflow from the Indian   (div <0)
364                  hdiv_160_89_kt(jk) = hdiv_160_89(jk) - zemp_red / (10. * e1v(ji,jj) * fse3v(ji,jj,jk) )
365               END DO
366               hdivn(ji, jj,:) = hdivn(ji, jj,:) + hdiv_160_89_kt(:)
367            END DO
368         END DO
369         !                  ! ---------------- !
370      CASE( 'tra' )         !  update (ta,sa)  ! (call by traadv module)
371         !                  ! --------=======- !
372         !
373         DO jj = mj0(88), mj1(88)              !** (161,88)   (Gulf of Aden side, north point)
374            DO ji = mi0(161), mi1(161) 
375               DO jk = 1, jpkm1                         ! surf inflow + reciculation (from Gulf of Aden)
376                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_161_88_kt(jk) * tn(ji,jj,jk)
377                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_161_88_kt(jk) * sn(ji,jj,jk)
378               END DO
379            END DO
380         END DO
381         DO jj = mj0(87), mj1(87)              !** (161,87)   (Gulf of Aden side, south point)
382            DO ji = mi0(161), mi1(161) 
383               jk =  21                                 ! deep outflow + recirulation (combined flux)
384               ta(ji,jj,jk) = ta(ji,jj,jk) + hdiv_161_88(20) * tn(ji  ,jj+1,20)   &  ! upper recirculation from Gulf of Aden
385                  &                        + hdiv_161_88(21) * tn(ji  ,jj+1,21)   &  ! deep  recirculation from Gulf of Aden
386                  &                        + hdiv_160_89(16) * tn(ji-1,jj+2,16)      ! deep inflow from Red sea
387               sa(ji,jj,jk) = sa(ji,jj,jk) + hdiv_161_88(20) * sn(ji  ,jj+1,20)   &
388                  &                        + hdiv_161_88(21) * sn(ji  ,jj+1,21)   &
389                  &                        + hdiv_160_89(16) * sn(ji-1,jj+2,16)   
390            END DO
391         END DO
392         DO jj = mj0(89), mj1(89)              !** (161,88)   (Red sea side)
393            DO ji = mi0(160), mi1(160)
394               DO jk = 1, 14                            ! surface inflow (from Gulf of Aden)
395                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_160_89_kt(jk) * tn(ji+1,jj-1,jk)
396                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_160_89_kt(jk) * sn(ji+1,jj-1,jk)
397               END DO
398               !                                        ! deep    outflow (from Red sea)
399               ta(ji,jj,16) = ta(ji,jj,16) - hdiv_160_89(jk) * tn(ji,jj,jk)
400               sa(ji,jj,16) = sa(ji,jj,16) - hdiv_160_89(jk) * sn(ji,jj,jk)
401            END DO
402         END DO
403         !
404         !                  ! ---------------- !
405      CASE( 'spg' )         !  update (ua,va)  ! (call by dynspg module)
406         !                  ! --------=======- !
407         ! at this stage, (ua,va) are the after velocity, not the tendancy
408         ! compute the velocity from the divergence at T-point
409         !
410         DO jj = mj0(88), mj1(88)              !** (160,88)   (Gulf of Aden side, north point)
411            DO ji = mi0(160), mi1(160)                   ! 160, not 161 as it is a U-point)
412               ua(ji,jj,:) = - hdiv_161_88_kt(:) / ( e1t(ji+1,jj) * e2t(ji+1,jj) * fse3t(ji+1,jj,:) )   &
413                  &                              * e2u(ji,jj) * fse3u(ji,jj,:)
414            END DO
415         END DO
416         DO jj = mj0(87), mj1(87)              !** (160,87)   (Gulf of Aden side, south point)
417            DO ji = mi0(160), mi1(160)                   ! 160, not 161 as it is a U-point)
418               ua(ji,jj,:) = - hdiv_161_87(:) / ( e1t(ji+1,jj) * e2t(ji+1,jj) * fse3t(ji+1,jj,:) )   &
419                  &                           * e2u(ji,jj) * fse3u(ji,jj,:)
420            END DO
421         END DO
422         DO jj = mj0(88), mj1(88)              !** profile of divergence at (160,89)   (Red sea side)
423            DO ji = mi0(160), mi1(160)                   ! 88, not 89 as it is a V-point)
424               va(ji,jj,:) = - hdiv_160_89_kt(:) / ( e1t(ji,jj+1) * e2t(ji,jj+1) * fse3t(ji,jj+1,:) )   &
425                  &                              * e1v(ji,jj) * fse3v(ji,jj,:)
426            END DO
427         END DO
428      END SELECT
429      !
430   END SUBROUTINE cla_bab_el_mandeb
431   
432
433   SUBROUTINE cla_gibraltar( cd_td )
434      !! -------------------------------------------------------------------
435      !!                 ***  ROUTINE cla_gibraltar  ***
436      !!       
437      !! ** Purpose :   update the now horizontal divergence, the tracer
438      !!              tendancyand the after velocity in vicinity of Gibraltar
439      !!              strait ( Persian Gulf - Indian ocean ).
440      !!
441      !! ** Method :
442      !!                     _______________________           
443      !!      deep  zio_flow /====|///////|====> surf. zio_flow
444      !!    + deep  zrecirc  \----|///////|     (+balance of emp)
445      !! 102                      u///////u
446      !!      mid.  recicul    <--|///////|<==== deep  zio_flow
447      !!                     _____|_______|_____ 
448      !!      surf. zio_flow ====>|///////|       
449      !!    (+balance of emp)     |///////|
450      !! 101                      u///////|             
451      !!      mid.  recicul    -->|///////|               Caution: zrecirc split into
452      !!      deep  zrecirc  ---->|///////|                  upper & bottom recirculation
453      !!                   _______|_______|_______
454      !!                     139     140     141 
455      !!
456      !!---------------------------------------------------------------------
457      CHARACTER(len=1), INTENT(in) ::   cd_td   ! ='div' update the divergence
458      !                                         ! ='tra' update the tracers
459      !                                         ! ='spg' update after velocity
460      INTEGER  ::   ji, jj, jk   ! dummy loop indices
461      REAL(wp) ::   zemp_med     ! temporary scalar
462      REAL(wp) ::   zio_flow, zrecirc_upp, zrecirc_mid, zrecirc_bot
463      !!---------------------------------------------------------------------
464      !
465      SELECT CASE( cd_td ) 
466      !                     ! ---------------- !
467      CASE( 'ini' )         !  initialisation  !
468         !                  ! ---------------- !
469         !                                     !** initialization of the velocity
470         hdiv_139_101(:) = 0.e0                     !  139,101 (Atlantic side, south point)
471         hdiv_139_102(:) = 0.e0                     !  139,102 (Atlantic side, north point)
472         hdiv_141_102(:) = 0.e0                     !  141,102 (Med sea  side)
473           
474         !                                     !** imposed transport
475         zio_flow    = 0.8e6                        ! inflow surface  water
476         zrecirc_mid = 0.7e6                        ! middle recirculation water
477         zrecirc_upp = 2.5e6                        ! upper  recirculation water
478         zrecirc_bot = 3.5e6                        ! bottom recirculation water
479         !
480         DO jj = mj0(101), mj1(101)            !** profile of hdiv at 139,101 (Atlantic side, south point)
481            DO ji = mi0(139), mi1(139)         !-----------------------------
482               DO jk = 1, 14                        ! surface in/out flow (Atl -> Med)   (div >0)
483                  hdiv_139_101(jk) = + zio_flow / ( 14. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
484               END DO
485               DO jk = 15, 20                       ! middle  reciculation (Atl 101 -> Atl 102)   (div >0)   
486                  hdiv_139_101(jk) = + zrecirc_mid / ( 6. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
487               END DO
488               !                                    ! upper reciculation (Atl 101 -> Atl 101)   (div >0)
489               hdiv_139_101(21) =               + zrecirc_upp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
490               !
491               !                                    ! upper & bottom reciculation (Atl 101 -> Atl 101 & 102)   (div >0)
492               hdiv_139_101(22) = ( zrecirc_bot - zrecirc_upp ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
493            END DO
494         END DO
495         DO jj = mj0(102), mj1(102)            !** profile of hdiv at 139,102 (Atlantic side, north point)
496            DO ji = mi0(139), mi1(139)         !-----------------------------
497               DO jk = 15, 20                       ! middle reciculation (Atl 101 -> Atl 102)   (div <0)               
498                  hdiv_139_102(jk) = - zrecirc_mid / ( 6. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
499               END DO
500               !                                    ! outflow of Mediterranean sea + deep recirculation   (div <0)
501               hdiv_139_102(22) = - ( zio_flow + zrecirc_bot ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
502            END DO
503         END DO
504         DO jj = mj0(102), mj1(102)            !** velocity profile at 141,102  (Med sea side)
505            DO ji = mi0(141), mi1(141)         !------------------------------
506               DO  jk = 1, 14                       ! surface inflow in the Med     (div <0)
507                  hdiv_141_102(jk) = - zio_flow / ( 14. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
508               END DO
509               !                                    ! deep    outflow toward the Atlantic    (div >0)
510               hdiv_141_102(21)    = + zio_flow / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
511            END DO
512         END DO
513         !                  ! ---------------- !
514      CASE( 'div' )         !   update hdivn   ! (call by divcur module)
515         !                  ! ---------=====-- !
516         !                                     !** emp on the Mediterranean Sea  (div >0)
517         zemp_med = 0.e0                       !-------------------------------
518         DO jj = mj0(96), mj1(110)                  ! sum over the Med sea
519            DO ji = mi0(141),mi1(181)
520               zemp_med = zemp_med + emp(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj) 
521            END DO
522         END DO
523         DO jj = mj0(96), mj1(96)                   ! minus 2 points in Red Sea
524            DO ji = mi0(148),mi1(148)
525               zemp_med = zemp_med - emp(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
526            END DO
527            DO ji = mi0(149),mi1(149)
528               zemp_med = zemp_med - emp(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
529            END DO
530         END DO
531         IF( lk_mpp )   CALL mpp_sum( zemp_med )    ! sum with other processors value
532         zemp_med = zemp_med * 1.e-3                ! convert in m3
533         !
534         !                                     !** Correct hdivn (including emp adjustment)
535         !                                     !-------------------------------------------
536         DO jj = mj0(101), mj1(101)                 !* 139,101 (Atlantic side, south point)
537            DO ji = mi0(139), mi1(139) 
538               hdiv_139_101_kt(:) = hdiv_139_101(:)     
539               DO jk = 1, 14                              ! increase the inflow from the Atlantic   (div >0)
540                  hdiv_139_101_kt(jk) = hdiv_139_101(jk) + zemp_med / ( 14. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
541               END DO
542               hdivn(ji, jj,:) = hdivn(ji, jj,:) + hdiv_139_101_kt(:)
543            END DO
544         END DO
545         DO jj = mj0(102), mj1(102)                 !* 139,102 (Atlantic side, north point)
546            DO ji = mi0(139), mi1(139) 
547               hdivn(ji,jj,:) = hdivn(ji,jj,:) + hdiv_139_102(:)
548            END DO
549         END DO
550         DO jj = mj0(102), mj1(102)                 !* 141,102 (Med side)
551            DO ji = mi0(141), mi1(141) 
552               hdiv_141_102(:) = hdiv_141_102(:)
553               DO jk = 1, 14                              ! increase the inflow from the Atlantic   (div <0)
554                  hdiv_141_102_kt(jk) = hdiv_141_102(jk) - zemp_med / ( 14. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
555               END DO
556               hdivn(ji, jj,:) = hdivn(ji, jj,:) + hdiv_141_102_kt(:)
557            END DO
558         END DO
559         !                  ! ---------------- !
560      CASE( 'tra' )         !  update (ta,sa)  ! (call by traadv module)
561         !                  ! --------=======- !
562         !
563         DO jj = mj0(101), mj1(101)            !** 139,101 (Atlantic side, south point)   (div >0)
564            DO ji = mi0(139), mi1(139) 
565               DO jk = 1, jpkm1                         ! surf inflow + mid. & bottom reciculation (from Atlantic)   
566                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_139_101_kt(jk) * tn(ji,jj,jk)
567                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_139_101_kt(jk) * sn(ji,jj,jk)
568               END DO
569            END DO
570         END DO
571         !
572         DO jj = mj0(102), mj1(102)            !** 139,102 (Atlantic side, north point)   (div <0)
573            DO ji = mi0(139), mi1(139) 
574               DO jk = 15, 20                            ! middle  reciculation (Atl 101 -> Atl 102)   (div <0)
575                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_139_102(jk) * tn(ji,jj-1,jk)  ! middle Atlantic recirculation
576                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_139_102(jk) * sn(ji,jj-1,jk)
577               END DO
578               !                                         ! upper & bottom Atl. reciculation (Atl 101 -> Atl 102) - (div <0)
579               !                                         ! deep Med flow                    (Med 102 -> Atl 102) - (div <0)
580               ta(ji,jj,22) = ta(ji,jj,22) + hdiv_141_102(21) * tn(ji+2,jj  ,21)   &  ! deep Med flow 
581                  &                        + hdiv_139_101(21) * tn(ji  ,jj-1,21)   &  ! upper  Atlantic recirculation 
582                  &                        + hdiv_139_101(22) * tn(ji  ,jj-1,22)      ! bottom Atlantic recirculation 
583               sa(ji,jj,22) = sa(ji,jj,22) + hdiv_141_102(21) * sn(ji+2,jj  ,21)   &
584                  &                        + hdiv_139_101(21) * sn(ji  ,jj-1,21)   &
585                  &                        + hdiv_139_101(22) * sn(ji  ,jj-1,22) 
586            END DO
587         END DO
588         DO jj = mj0(102), mj1(102)                 !* 141,102 (Med side)   (div <0)
589            DO ji = mi0(141), mi1(141) 
590               DO jk = 1, 14                             ! surface flow from Atlantic to Med sea
591                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_141_102_kt(jk) * tn(ji-2,jj-1,jk)
592                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_141_102_kt(jk) * sn(ji-2,jj-1,jk)
593               END DO
594               !                                         ! deeper flow from Med sea to Atlantic
595               ta(ji,jj,21) = ta(ji,jj,21) - hdiv_141_102(21) * tn(ji,jj,21)
596               sa(ji,jj,21) = sa(ji,jj,21) - hdiv_141_102(21) * sn(ji,jj,21)
597            END DO
598         END DO
599         !                  ! ---------------- !
600      CASE( 'spg' )         !  update (ua,va)  ! (call by dynspg module)
601         !                  ! --------=======- !
602         ! at this stage, (ua,va) are the after velocity, not the tendancy
603         ! compute the velocity from the divergence at T-point
604         !
605         DO jj = mj0(101), mj1(101)            !** 139,101 (Atlantic side, south point)
606            DO ji = mi0(139), mi1(139)                    ! div >0 => ua >0, same sign
607               ua(ji,jj,:) = hdiv_139_101_kt(:) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,:) )   &
608                  &                             * e2u(ji,jj) * fse3u(ji,jj,:)
609            END DO
610         END DO
611         DO jj = mj0(102), mj1(102)            !** 139,102 (Atlantic side, north point)
612            DO ji = mi0(139), mi1(139)                    ! div <0 => ua <0, same sign
613               ua(ji,jj,:) = hdiv_139_102(:) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,:) )   &
614                  &                          * e2u(ji,jj) * fse3u(ji,jj,:)   
615            END DO
616         END DO
617         DO jj = mj0(102), mj1(102)            !** 140,102 (Med side) (140 not 141 as it is a U-point)
618            DO ji = mi0(140), mi1(140)                    ! div >0 => ua <0, opposite sign
619               ua(ji,jj,:) = - hdiv_141_102(:) / ( e1t(ji+1,jj) * e2t(ji+1,jj) * fse3t(ji+1,jj,:) )   &
620                  &                            * e2u(ji,jj) * fse3u(ji,jj,:)
621            END DO
622         END DO
623         !
624      END SELECT
625      !
626   END SUBROUTINE cla_gibraltar
627
628
629   SUBROUTINE cla_hormuz( cd_td )
630      !! -------------------------------------------------------------------
631      !!                   ***  ROUTINE div_hormuz  ***
632      !!             
633      !! ** Purpose :   update the now horizontal divergence, the tracer
634      !!              tendancyand the after velocity in vicinity of Hormuz
635      !!              strait ( Persian Gulf - Indian ocean ).
636      !!
637      !! ** Method  :   Hormuz strait
638      !!            ______________   
639      !!            |/////|<==      surface inflow
640      !!        94  |/////|     
641      !!            |/////|==>      deep    outflow
642      !!            |_____|_______
643      !!              171    172     
644      !!---------------------------------------------------------------------
645      CHARACTER(len=1), INTENT(in) ::   cd_td   ! ='ini' initialisation
646      !!                                        ! ='div' update the divergence
647      !!                                        ! ='tra' update the tracers
648      !!                                        ! ='spg' update after velocity
649      !!
650      INTEGER  ::   ji, jj, jk   ! dummy loop indices
651      REAL(wp) ::   zio_flow     ! temporary scalar
652      !!---------------------------------------------------------------------
653      !
654      SELECT CASE( cd_td ) 
655      !                     ! ---------------- !
656      CASE( 'ini' )         !  initialisation  !
657         !                  ! ---------------- !
658         !                                     !** profile of horizontal divergence due to cross-land advection
659         zio_flow  = 1.e6                          ! imposed in/out flow
660         !
661         hdiv_172_94(:) = 0.e0         
662         !
663         DO jj = mj0(94), mj1(94)                  ! in/out flow at (i,j) = (172,94)
664            DO ji = mi0(172), mi1(172) 
665               DO jk = 1, 8                            ! surface inflow  (Indian ocean to Persian Gulf) (div<0)
666                  hdiv_172_94(jk) = - ( zio_flow / 8.e0 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
667               END DO
668               DO jk = 16, 18                          ! deep    outflow (Persian Gulf to Indian ocean) (div>0)
669                  hdiv_172_94(jk) = + ( zio_flow / 3.e0 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
670               END DO
671            END DO
672         END DO
673         !                                     !** T & S profile in the Hormuz strait (use in deep outflow)
674         !      Temperature       and         Salinity
675         t_171_94_hor(:)  = 0.e0   ;   s_171_94_hor(:)  = 0.e0
676         t_171_94_hor(16) = 18.4   ;   s_171_94_hor(16) = 36.27
677         t_171_94_hor(17) = 17.8   ;   s_171_94_hor(17) = 36.4
678         t_171_94_hor(18) = 16.    ;   s_171_94_hor(18) = 36.27
679         !
680         !                  ! ---------------- !
681      CASE( 'div' )         !   update hdivn   ! (call by divcur module)
682         !                  ! ---------=====-- !
683         !                                   
684         DO jj = mj0(94), mj1(94)              !** 172,94 (Indian ocean side)
685            DO ji = mi0(172), mi1(172) 
686               hdivn(ji,jj,:) = hdivn(ji,jj,:) + hdiv_172_94(:)
687            END DO
688         END DO
689         !                  ! ---------------- !
690      CASE( 'tra' )         !  update (ta,sa)  ! (call by traadv module)
691         !                  ! --------=======- !
692         !                         
693         DO jj = mj0(94), mj1(94)              !** 172,94 (Indian ocean side)
694            DO ji = mi0(172), mi1(172) 
695               DO jk = 1, 8                          ! surface inflow   (Indian ocean to Persian Gulf) (div<0)
696                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_172_94(jk) * tn(ji,jj,jk) 
697                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_172_94(jk) * sn(ji,jj,jk) 
698               END DO
699               DO jk = 16, 18                        ! deep outflow     (Persian Gulf to Indian ocean) (div>0)
700                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_172_94(jk) * t_171_94_hor(jk)
701                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_172_94(jk) * s_171_94_hor(jk)
702               END DO
703            END DO
704         END DO
705         !                  ! ---------------- !
706      CASE( 'spg' )         !  update (ua,va)  ! (call by dynspg module)
707         !                  ! --------=======- !
708         ! No barotropic flow through Hormuz strait
709         ! at this stage, (ua,va) are the after velocity, not the tendancy
710         ! compute the velocity from the divergence at T-point
711         DO jj = mj0(94), mj1(94)              !** 171,94 (Indian ocean side) (171 not 172 as it is the western U-point)
712            DO ji = mi0(171), mi1(171)                ! div >0 => ua >0, opposite sign
713               ua(ji,jj,:) = - hdiv_172_94(:) / ( e1t(ji+1,jj) * e2t(ji+1,jj) * fse3t(ji+1,jj,:) )   &
714                  &                           * e2u(ji,jj) * fse3u(ji,jj,:)
715            END DO
716         END DO
717         !
718      END SELECT
719      !
720   END SUBROUTINE cla_hormuz
721   
722#else
723   !!----------------------------------------------------------------------
724   !!   Default key                                            Dummy module
725   !!----------------------------------------------------------------------
726   USE lib_mpp, ONLY:   ctl_stop
727CONTAINS
728   SUBROUTINE cla_init
729      CALL ctl_stop( 'cla_init: Cross Land Advection hard coded for ORCA_R2 with 31 levels' )
730   END SUBROUTINE cla_init
731   SUBROUTINE cla_div( kt )
732      WRITE(*,*) 'cla_div: You should have not see this print! error?', kt
733   END SUBROUTINE cla_div
734   SUBROUTINE cla_traadv( kt ) 
735      WRITE(*,*) 'cla_traadv: You should have not see this print! error?', kt
736   END SUBROUTINE cla_traadv
737   SUBROUTINE cla_dynspg( kt ) 
738      WRITE(*,*) 'dyn_spg_cla: You should have not see this print! error?', kt
739   END SUBROUTINE cla_dynspg
740#endif
741   
742   !!======================================================================
743END MODULE cla
Note: See TracBrowser for help on using the repository browser.