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/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/LBC – NEMO

source: branches/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/LBC/cla.F90 @ 10207

Last change on this file since 10207 was 10207, checked in by cmao, 6 years ago

remove svn keyword

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