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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/cla.F90 @ 4429

Last change on this file since 4429 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

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