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.
Changeset 2392 for branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla.F90 – NEMO

Ignore:
Timestamp:
2010-11-15T22:20:05+01:00 (13 years ago)
Author:
gm
Message:

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla.F90

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