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 1438 for trunk/NEMO/OPA_SRC/DOM/domvvl.F90 – NEMO

Ignore:
Timestamp:
2009-05-11T16:34:47+02:00 (15 years ago)
Author:
rblod
Message:

Merge VVL branch with the trunk (act II), see ticket #429

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/DOM/domvvl.F90

    r1146 r1438  
    44   !! Ocean :  
    55   !!====================================================================== 
    6    !! History :   9.0  !  06-06  (B. Levier, L. Marie)  original code 
    7    !!              "   !  07-07  (D. Storkey) Bug fixes and code for BDY option. 
     6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code 
     7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate 
    88   !!---------------------------------------------------------------------- 
    9  
     9#if defined key_vvl 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   'key_vvl'                              variable volume 
    1212   !!---------------------------------------------------------------------- 
     13   !!   dom_vvl     : defined coefficients to distribute ssh on each layers 
    1314   !!---------------------------------------------------------------------- 
    14    !!   dom_vvl         : defined scale factors & depths at each time step 
    15    !!   dom_vvl_ini     : defined coefficients to distribute ssh on each layers 
    16    !!   dom_vvl_ssh     : defined the ocean sea level at each time step 
    17    !!---------------------------------------------------------------------- 
    18    !! * Modules used 
    1915   USE oce             ! ocean dynamics and tracers 
    2016   USE dom_oce         ! ocean space and time domain 
    2117   USE sbc_oce         ! surface boundary condition: ocean 
    22    USE dynspg_oce      ! surface pressure gradient variables 
    2318   USE phycst          ! physical constants 
    2419   USE in_out_manager  ! I/O manager 
    2520   USE lib_mpp         ! distributed memory computing library 
    2621   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    27    USE bdy_oce         ! unstructured open boundary conditions 
    2822 
    2923   IMPLICIT NONE 
    3024   PRIVATE 
    3125 
    32    !! * Routine accessibility 
    33    PUBLIC dom_vvl_ini    ! called by dom_init.F90 
    34    PUBLIC dom_vvl_ssh    ! called by trazdf.F90 
    35    PUBLIC dom_vvl        ! called by istate.F90 and step.F90 
    36    PUBLIC dom_vvl_sf_ini !  
    37    PUBLIC dom_vvl_sf     !  
     26   PUBLIC dom_vvl        ! called by domain.F90 
    3827 
    39    !! * Module variables 
    40    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !: 
    41       mut, muu, muv, muf                            !: 
     28   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ee_t, ee_u, ee_v, ee_f   !: ??? 
     29 
     30   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   mut, muu, muv, muf   !: ???  
    4231 
    4332   REAL(wp), DIMENSION(jpk) ::   r2dt               ! vertical profile time-step, = 2 rdttra  
     
    4837#  include "vectopt_loop_substitute.h90" 
    4938   !!---------------------------------------------------------------------- 
    50    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     39   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    5140   !! $Id$ 
    5241   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    5544CONTAINS        
    5645 
    57 #if defined key_vvl 
    58  
    59    SUBROUTINE dom_vvl_ini 
     46   SUBROUTINE dom_vvl 
    6047      !!---------------------------------------------------------------------- 
    61       !!                ***  ROUTINE dom_vvl_ini  *** 
     48      !!                ***  ROUTINE dom_vvl  *** 
    6249      !!                    
    6350      !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread 
     
    6552      !! 
    6653      !!---------------------------------------------------------------------- 
    67       INTEGER  ::   ji, jj, jk, zpk 
    68       REAL(wp), DIMENSION(jpi,jpj) ::   zmut, zmuu, zmuv, zmuf   ! 2D workspace 
     54      INTEGER  ::   ji, jj, jk 
     55      REAL(wp) ::   zcoefu, zcoefv, zcoeff 
    6956      !!---------------------------------------------------------------------- 
    7057 
    7158      IF(lwp)   THEN 
    7259         WRITE(numout,*) 
    73          WRITE(numout,*) 'dom_vvl_ini : Variable volume activated' 
    74          WRITE(numout,*) '~~~~~~~~~~~   compute coef. used to spread ssh over each layers' 
     60         WRITE(numout,*) 'dom_vvl : Variable volume activated' 
     61         WRITE(numout,*) '~~~~~~~~  compute coef. used to spread ssh over each layers' 
    7562      ENDIF 
    76  
    77       IF( ln_zps )  CALL ctl_stop( 'dom_vvl_ini : option ln_zps is incompatible with variable volume option key_vvl') 
    7863 
    7964#if defined key_zco  ||  defined key_dynspg_rl 
     
    8166#endif 
    8267 
    83       fsvdept (:,:,:) = gdept (:,:,:) 
    84       fsvdepw (:,:,:) = gdepw (:,:,:) 
    85       fsvde3w (:,:,:) = gdep3w(:,:,:) 
    86       fsve3t (:,:,:) = e3t   (:,:,:) 
    87       fsve3u (:,:,:) = e3u   (:,:,:) 
    88       fsve3v (:,:,:) = e3v   (:,:,:) 
    89       fsve3f (:,:,:) = e3f   (:,:,:) 
    90       fsve3w (:,:,:) = e3w   (:,:,:) 
    91       fsve3uw (:,:,:) = e3uw  (:,:,:) 
    92       fsve3vw (:,:,:) = e3vw  (:,:,:) 
     68      fsdept(:,:,:) = gdept (:,:,:) 
     69      fsdepw(:,:,:) = gdepw (:,:,:) 
     70      fsde3w(:,:,:) = gdep3w(:,:,:) 
     71      fse3t (:,:,:) = e3t   (:,:,:) 
     72      fse3u (:,:,:) = e3u   (:,:,:) 
     73      fse3v (:,:,:) = e3v   (:,:,:) 
     74      fse3f (:,:,:) = e3f   (:,:,:) 
     75      fse3w (:,:,:) = e3w   (:,:,:) 
     76      fse3uw(:,:,:) = e3uw  (:,:,:) 
     77      fse3vw(:,:,:) = e3vw  (:,:,:) 
    9378 
    9479      ! mu computation 
    95       zmut(:,:)   = 0.e0 
    96       zmuu(:,:)   = 0.e0 
    97       zmuv(:,:)   = 0.e0 
    98       zmuf(:,:)   = 0.e0 
    99       mut (:,:,:) = 0.e0 
    100       muu (:,:,:) = 0.e0 
    101       muv (:,:,:) = 0.e0 
    102       muf (:,:,:) = 0.e0 
     80      ! -------------- 
     81      ! define ee_t, u, v and f as in sigma coordinate (ee_t = 1/ht, ...) 
     82      ee_t(:,:) = fse3t_0(:,:,1)        ! Lower bound : thickness of the first model level 
     83      ee_u(:,:) = fse3u_0(:,:,1) 
     84      ee_v(:,:) = fse3v_0(:,:,1) 
     85      ee_f(:,:) = fse3f_0(:,:,1) 
     86      DO jk = 2, jpkm1                   ! Sum of the masked vertical scale factors 
     87         ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 
     88         ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
     89         ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
     90         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask 
     91            ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk) 
     92         END DO 
     93      END DO   
     94      !                                  ! Compute and mask the inverse of the local depth at T, U, V and F points 
     95      ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1) 
     96      ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1) 
     97      ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1) 
     98      DO jj = 1, jpjm1                         ! f-point case fmask cannot be used  
     99         ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 
     100      END DO 
     101      CALL lbc_lnk( ee_f, 'F', 1. )       ! lateral boundary condition on ee_f 
     102      ! 
     103      DO jk = 1, jpk 
     104         mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk)   ! at T levels 
     105         muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk)   ! at T levels 
     106         muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk)   ! at T levels 
     107      END DO 
     108      DO jk = 1, jpk 
     109         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask 
     110               muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels 
     111         END DO 
     112         muf(:,jpj,jk) = 0.e0 
     113      END DO 
     114      CALL lbc_lnk( muf, 'F', 1. )       ! lateral boundary condition on ee_f 
    103115 
    104       DO jj = 1, jpj 
    105          DO ji = 1, jpi 
    106             zpk = mbathy(ji,jj) - 1 
    107             DO jk = 1, zpk 
    108                zmut(ji,jj) = zmut(ji,jj) + fsve3t(ji,jj,jk) * SUM( fsve3t(ji,jj,jk:zpk) ) 
    109                zmuu(ji,jj) = zmuu(ji,jj) + fsve3u(ji,jj,jk) * SUM( fsve3u(ji,jj,jk:zpk) ) 
    110                zmuv(ji,jj) = zmuv(ji,jj) + fsve3v(ji,jj,jk) * SUM( fsve3v(ji,jj,jk:zpk) ) 
    111                zmuf(ji,jj) = zmuf(ji,jj) + fsve3f(ji,jj,jk) * SUM( fsve3f(ji,jj,jk:zpk) ) 
    112             END DO 
     116 
     117      ! Reference ocean depth at U- and V-points 
     118      hu_0(:,:) = 0.e0     
     119      hv_0(:,:) = 0.e0 
     120      DO jk = 1, jpk 
     121         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
     122         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
     123      END DO 
     124 
     125      ! before and now Sea Surface Height at u-, v-, f-points 
     126      DO jj = 1, jpjm1 
     127         DO ji = 1, jpim1 
     128            zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     129            zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     130            zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
     131            sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
     132               &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )    
     133            sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
     134               &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )    
     135            sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 & 
     136               &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )                
     137            sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
     138               &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )    
     139            sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &  
     140               &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )      
     141            sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 & 
     142               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )                
    113143         END DO 
    114144      END DO 
    115  
    116       DO jj = 1, jpj 
    117          DO ji = 1, jpi 
    118             zpk = mbathy(ji,jj) - 1 
    119             DO jk = 1, zpk 
    120 #if defined key_sigma_vvl 
    121                mut(ji,jj,jk) = 1./SUM( fsve3t(ji,jj,1:zpk) )  
    122                muu(ji,jj,jk) = 1./SUM( fsve3u(ji,jj,1:zpk) )  
    123                muv(ji,jj,jk) = 1./SUM( fsve3v(ji,jj,1:zpk) )  
    124                muf(ji,jj,jk) = 1./SUM( fsve3f(ji,jj,1:zpk) )  
    125 #else 
    126                mut(ji,jj,jk) = SUM( fsve3t(ji,jj,jk:zpk) ) / zmut(ji,jj) 
    127                muu(ji,jj,jk) = SUM( fsve3u(ji,jj,jk:zpk) ) / zmuu(ji,jj) 
    128                muv(ji,jj,jk) = SUM( fsve3v(ji,jj,jk:zpk) ) / zmuv(ji,jj) 
    129                muf(ji,jj,jk) = SUM( fsve3f(ji,jj,jk:zpk) ) / zmuf(ji,jj) 
    130 #endif 
    131             END DO 
    132          END DO 
    133       END DO 
    134  
    135  
    136    END SUBROUTINE dom_vvl_ini 
    137  
    138  
    139  
    140    SUBROUTINE dom_vvl 
    141       !!---------------------------------------------------------------------- 
    142       !!                ***  ROUTINE dom_vvl  *** 
    143       !!                    
    144       !! ** Purpose :  compute ssh at U-V-F points, T-W scale factors and local 
    145       !!               depths at each time step. 
    146       !!---------------------------------------------------------------------- 
    147       !! * Local declarations 
    148       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    149       REAL(wp), DIMENSION(jpi,jpj) :: zsshf    ! 2D workspace 
    150       !!---------------------------------------------------------------------- 
    151  
    152       ! Sea Surface Height at u-v-fpoints 
    153       DO jj = 1, jpjm1 
    154          DO ji = 1,jpim1 
    155             sshu(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  & 
    156                &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    157                &            + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
    158             ! 
    159             sshv(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  & 
    160                &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
    161                &            + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
    162             ! 
    163             zsshf(ji,jj) = 0.25 * fmask(ji,jj,1)                 & 
    164                &           * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
    165                &             + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 
    166          END DO 
    167       END DO 
    168  
    169145      ! Boundaries conditions 
    170       CALL lbc_lnk( sshu,  'U', 1. ) 
    171       CALL lbc_lnk( sshv,  'V', 1. ) 
    172       CALL lbc_lnk( zsshf, 'F', 1. ) 
    173  
    174       ! Scale factors at T levels 
    175       DO jk = 1, jpkm1 
    176          fse3t(:,:,jk) = fsve3t(:,:,jk) * ( 1 + sshn(:,:)  * mut(:,:,jk) ) 
    177          fse3u(:,:,jk) = fsve3u(:,:,jk) * ( 1 + sshu(:,:)  * muu(:,:,jk) ) 
    178          fse3v(:,:,jk) = fsve3v(:,:,jk) * ( 1 + sshv(:,:)  * muv(:,:,jk) ) 
    179          fse3f(:,:,jk) = fsve3f(:,:,jk) * ( 1 + zsshf(:,:) * muf(:,:,jk) ) 
    180       END DO 
    181  
    182       ! Scale factors at W levels 
    183       fse3w (:,:,1) = fse3t(:,:,1) 
    184       fse3uw(:,:,1) = fse3u(:,:,1) 
    185       fse3vw(:,:,1) = fse3v(:,:,1) 
    186       DO jk = 2, jpk 
    187          fse3w (:,:,jk) = 0.5 * ( fse3t(:,:,jk-1) + fse3t(:,:,jk) ) 
    188          fse3uw(:,:,jk) = 0.5 * ( fse3u(:,:,jk-1) + fse3u(:,:,jk) ) 
    189          fse3vw(:,:,jk) = 0.5 * ( fse3v(:,:,jk-1) + fse3v(:,:,jk) ) 
    190       END DO 
    191  
    192       ! T and W points depth 
    193       fsdept(:,:,1) = 0.5 * fse3w(:,:,1) 
    194       fsdepw(:,:,1) = 0.e0 
    195       fsde3w(:,:,1) = fsdept(:,:,1) - sshn(:,:) 
    196       DO jk = 2, jpk 
    197          fsdept(:,:,jk) = fsdept(:,:,jk-1) + fse3w(:,:,jk) 
    198          fsdepw(:,:,jk) = fsdepw(:,:,jk-1) + fse3t(:,:,jk-1) 
    199          fsde3w(:,:,jk) = fsdept(:,:,jk  ) - sshn(:,:) 
    200       END DO 
    201  
    202       IF( MINVAL(fse3t(:,:,:)) < 0.0 ) THEN 
    203          CALL ctl_stop('E R R O R : dom_vvl','Level thickness fse3t less than zero.') 
    204          nstop = nstop+1 
    205       ENDIF 
    206  
    207       ! Local depth or Inverse of the local depth of the water column at u- and v-points 
    208       ! ------------------------------ 
    209  
    210       ! Ocean depth at U- and V-points 
    211       hu(:,:) = 0.e0    ;    hv(:,:) = 0.e0 
    212  
    213       DO jk = 1, jpk 
    214          hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    215          hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
    216       END DO 
    217  
    218  
    219       ! Inverse of the local depth 
    220       hur(:,:) = fse3u(:,:,1)             ! Lower bound : thickness of the first model level 
    221       hvr(:,:) = fse3v(:,:,1) 
    222        
    223       DO jk = 2, jpk                      ! Sum of the vertical scale factors 
    224          hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    225          hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
    226       END DO 
    227  
    228       ! Compute and mask the inverse of the local depth 
    229       hur(:,:) = 1. / hur(:,:) * umask(:,:,1) 
    230       hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1) 
    231  
     146      CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
     147      CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
     148      CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. ) 
     149      ! 
    232150   END SUBROUTINE dom_vvl 
    233151 
    234  
    235  
    236    SUBROUTINE dom_vvl_ssh( kt )  
    237       !!---------------------------------------------------------------------- 
    238       !!                ***  ROUTINE dom_vvl_ssh  *** 
    239       !!                    
    240       !! ** Purpose :  compute the ssha field used in tra_zdf, dynnxt, tranxt  
    241       !!               and dynspg routines 
    242       !!---------------------------------------------------------------------- 
    243       !! * Arguments 
    244       INTEGER, INTENT( in )  ::    kt                         ! time step 
    245  
    246       !! * Local declarations 
    247       INTEGER  ::   ji, jj, jk                                ! dummy loop indices 
    248       REAL(wp) ::   z2dt, zraur                               ! temporary scalars 
    249       REAL(wp), DIMENSION(jpi,jpj) ::   zun, zvn, zhdiv       ! 2D workspace 
    250       !!---------------------------------------------------------------------- 
    251  
    252       !! Sea surface height after (ssha) in variable volume case 
    253       !! --------------------------====-----===============----- 
    254       !! ssha is needed for the tracer time-stepping (trazdf_imp or 
    255       !! tra_nxt), for the ssh time-stepping (dynspg_*) and for the 
    256       !! dynamics time-stepping (dynspg_flt or dyn_nxt, and wzv). 
    257       !! un, vn (or un_b and vn_b) and emp are needed, so it must be 
    258       !! done after the call of oce_sbc. 
    259  
    260       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000 
    261          r2dt(:) =  rdttra(:)                       ! = rdtra (restarting with Euler time stepping) 
    262       ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1 
    263          r2dt(:) = 2. * rdttra(:)                   ! = 2 rdttra (leapfrog) 
    264       ENDIF 
    265  
    266       z2dt  = r2dt(1) 
    267       zraur = 1. / rauw 
    268  
    269       ! Vertically integrated quantities 
    270       ! -------------------------------- 
    271       IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN 
    272          zun(:,:) = 0.e0    ;    zvn(:,:) = 0.e0 
    273          ! 
    274          DO jk = 1, jpkm1 
    275             !                                               ! Vertically integrated transports (now) 
    276             zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    277             zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    278          END DO 
    279       ELSE ! lk_dynspg_ts is true 
    280          zun (:,:) = un_b(:,:)    ;    zvn (:,:) = vn_b(:,:) 
    281       ENDIF 
    282  
    283       ! Horizontal divergence of barotropic transports 
    284       !-------------------------------------------------- 
    285       DO jj = 2, jpjm1 
    286          DO ji = 2, jpim1   ! vector opt. 
    287             zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )    & 
    288                &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )    & 
    289                &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )    & 
    290                &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )  & 
    291                &           / ( e1t(ji,jj) * e2t(ji,jj) ) 
    292          END DO 
    293       END DO 
    294  
    295 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 
    296       ! open boundaries (div must be zero behind the open boundary) 
    297       !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    298       IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east 
    299       IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west 
    300       IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north 
    301       IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south 
    302 #endif 
    303  
    304 #if defined key_bdy 
    305       zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) 
    306 #endif 
    307  
    308       CALL lbc_lnk( zhdiv, 'T', 1. ) 
    309  
    310       ! Sea surface elevation time stepping 
    311       ! ----------------------------------- 
    312       IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN 
    313          ssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 
    314       ELSE ! lk_dynspg_ts is true 
    315          ssha(:,:) = (  sshb_b(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    316       ENDIF 
    317  
    318  
    319    END SUBROUTINE dom_vvl_ssh 
    320  
    321    SUBROUTINE  dom_vvl_sf( zssh, gridp, sfe3 ) 
    322       !!---------------------------------------------------------------------- 
    323       !!                ***  ROUTINE dom_vvl_sf  *** 
    324       !!                    
    325       !! ** Purpose :  compute vertical scale factor at each time step 
    326       !!---------------------------------------------------------------------- 
    327       !! * Arguments 
    328       CHARACTER(LEN=1)                , INTENT( in ) :: gridp   ! grid point type 
    329       REAL(wp), DIMENSION(jpi,jpj)    , INTENT( in ) :: zssh    ! 2D workspace 
    330       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: sfe3    ! 3D workspace 
    331  
    332       !! * Local declarations 
    333       INTEGER  ::   jk                                      ! dummy loop indices 
    334       !!---------------------------------------------------------------------- 
    335  
    336       SELECT CASE (gridp) 
    337  
    338       CASE ('U') 
    339          ! 
    340          DO jk = 1, jpk 
    341             sfe3(:,:,jk)  = fsve3u(:,:,jk) * ( 1 + zssh(:,:) * muu(:,:,jk) ) 
    342          ENDDO 
    343  
    344       CASE ('V') 
    345          ! 
    346          DO jk = 1, jpk 
    347             sfe3(:,:,jk)  = fsve3v(:,:,jk) * ( 1 + zssh(:,:) * muv(:,:,jk) ) 
    348          ENDDO 
    349  
    350       CASE ('T') 
    351          ! 
    352          DO jk = 1, jpk 
    353             sfe3(:,:,jk)  = fsve3t(:,:,jk) * ( 1 + zssh(:,:) * mut(:,:,jk) ) 
    354          ENDDO 
    355  
    356       END SELECT 
    357  
    358    END SUBROUTINE dom_vvl_sf 
    359  
    360    SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini ) 
    361       !!---------------------------------------------------------------------- 
    362       !!                ***  ROUTINE dom_vvl_sf_ini  *** 
    363       !!                    
    364       !! ** Purpose :  affect the appropriate vertical scale factor. It is done 
    365       !!               to avoid compiling error when using key_zco cpp key 
    366       !!---------------------------------------------------------------------- 
    367       !! * Arguments 
    368       CHARACTER(LEN=1)                , INTENT( in ) :: gridp      ! grid point type 
    369       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini    ! 3D workspace 
    370       !!---------------------------------------------------------------------- 
    371  
    372       SELECT CASE (gridp) 
    373  
    374       CASE ('U') 
    375          ! 
    376          sfe3ini(:,:,:)  = fse3u(:,:,:) 
    377  
    378       CASE ('V') 
    379          ! 
    380          sfe3ini(:,:,:)  = fse3v(:,:,:) 
    381  
    382       CASE ('T') 
    383          ! 
    384          sfe3ini(:,:,:)  = fse3t(:,:,:) 
    385  
    386       END SELECT 
    387  
    388    END SUBROUTINE dom_vvl_sf_ini 
    389152#else 
    390153   !!---------------------------------------------------------------------- 
    391154   !!   Default option :                                      Empty routine 
    392155   !!---------------------------------------------------------------------- 
    393    SUBROUTINE dom_vvl_ini 
    394    END SUBROUTINE dom_vvl_ini 
     156CONTAINS 
    395157   SUBROUTINE dom_vvl 
    396158   END SUBROUTINE dom_vvl 
    397    SUBROUTINE dom_vvl_ssh( kt )  
    398       !! * Arguments 
    399       INTEGER, INTENT( in )  ::    kt        ! time step 
    400       WRITE(*,*) 'dom_vvl_ssh: You should not have seen this print! error?', kt 
    401    END SUBROUTINE dom_vvl_ssh 
    402    SUBROUTINE dom_vvl_sf( zssh, gridp, sfe3 )  
    403       !! * Arguments 
    404       CHARACTER(LEN=1)                , INTENT( in ) :: gridp   ! grid point type 
    405       REAL(wp), DIMENSION(jpi,jpj)    , INTENT( in ) :: zssh    ! 2D workspace 
    406       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3    ! 3D workspace 
    407       sfe3(:,:,:) = 0.e0 
    408       WRITE(*,*) 'sfe3: You should not have seen this print! error?', gridp 
    409       WRITE(*,*) 'sfe3: You should not have seen this print! error?', zssh(1,1) 
    410    END SUBROUTINE dom_vvl_sf 
    411    SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini )  
    412       !! * Arguments 
    413       CHARACTER(LEN=1)                , INTENT( in ) :: gridp    ! grid point type 
    414       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini  ! 3D workspace 
    415       sfe3ini(:,:,:) = 0.e0 
    416       WRITE(*,*) 'sfe3ini: You should not have seen this print! error?', gridp 
    417    END SUBROUTINE dom_vvl_sf_ini 
    418159#endif 
    419160 
Note: See TracChangeset for help on using the changeset viewer.