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 2370 for branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90 – NEMO

Ignore:
Timestamp:
2010-11-10T08:48:54+01:00 (13 years ago)
Author:
gm
Message:

v3.3beta: ice-ocean stress at kt with VP & EVP (LIM-2 and -3)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90

    r2332 r2370  
    44   !!   Sea-Ice dynamics :   
    55   !!====================================================================== 
    6    !! History :   1.0  !  01-04  (LIM)  Original code 
    7    !!             2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp 
    8    !!             2.0  !  03-08  (C. Ethe) add lim_dyn_init 
    9    !!             2.0  !  06-07  (G. Madec)  Surface module 
     6   !! History :  1.0  ! 2001-04  (LIM)  Original code 
     7   !!            2.0  ! 2002-08  (C. Ethe, G. Madec)  F90, mpp 
     8   !!            2.0  ! 2003-08  (C. Ethe) add lim_dyn_init 
     9   !!            2.0  ! 2006-07  (G. Madec)  Surface module 
     10   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case 
    1011   !!--------------------------------------------------------------------- 
    1112#if defined key_lim2 
     
    1617   !!    lim_dyn_init_2 : initialization and namelist read 
    1718   !!---------------------------------------------------------------------- 
    18    USE dom_oce        ! ocean space and time domain 
    19    USE sbc_oce        ! 
    20    USE phycst         ! 
    21    USE ice_2          ! 
    22    USE dom_ice_2      ! 
    23    USE limistate_2    ! 
    24 #if defined key_lim2_vp 
    25    USE limrhg_2       ! ice rheology 
    26 #else 
    27    USE limrhg         ! LIM : EVP ice rheology 
    28 #endif 
    29    USE lbclnk         ! 
    30    USE lib_mpp        ! 
    31    USE in_out_manager ! I/O manager 
    32    USE prtctl         ! Print control 
     19   USE dom_oce          ! ocean space and time domain 
     20   USE sbc_oce          ! ocean surface boundary condition 
     21   USE phycst           ! physical constant 
     22   USE ice_2            ! LIM-2: ice variables 
     23   USE sbc_ice          ! Surface boundary condition: sea-ice fields 
     24   USE dom_ice_2        ! LIM-2: ice domain 
     25   USE limistate_2      ! LIM-2: initial state 
     26   USE limrhg_2         ! LIM-2: VP  ice rheology 
     27   USE limrhg           ! LIM  : EVP ice rheology 
     28   USE lbclnk           ! lateral boundary condition - MPP link 
     29   USE lib_mpp          ! MPP library 
     30   USE in_out_manager   ! I/O manager 
     31   USE prtctl           ! Print control 
    3332 
    3433   IMPLICIT NONE 
    3534   PRIVATE 
    3635 
    37    PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim 
    38  
    39    !! * Module variables 
    40    REAL(wp)  ::  rone    = 1.e0   ! constant value 
    41  
     36   PUBLIC   lim_dyn_2   ! routine called by sbc_ice_lim 
     37 
     38   !! * Substitutions 
    4239#  include "vectopt_loop_substitute.h90" 
    4340   !!---------------------------------------------------------------------- 
    4441   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010) 
    4542   !! $Id$ 
    46    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    47    !!---------------------------------------------------------------------- 
    48  
     43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     44   !!---------------------------------------------------------------------- 
    4945CONTAINS 
    5046 
     
    9086            i_jpj = jpj 
    9187            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    92 #if defined key_lim2_vp 
    93             CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology 
    94 #else 
    95             CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology 
    96 #endif 
     88            IF( lk_lim2_vp )   THEN   ;   CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology 
     89            ELSE                      ;   CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology 
     90            ENDIF 
    9791            ! 
    9892         ELSE                                 ! optimization of the computational area 
     
    112106                  i_j1 = i_j1 + 1 
    113107               END DO 
    114 #if defined key_lim2_vp 
    115                i_j1 = MAX( 1, i_j1-1 ) 
    116                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    117                !  
    118                CALL lim_rhg_2( i_j1, i_jpj ) 
    119                !  
    120 #else 
    121                i_j1 = MAX( 1, i_j1-2 ) 
     108               IF( lk_lim2_vp )   THEN             ! VP  rheology 
     109                  i_j1 = MAX( 1, i_j1-1 ) 
     110                  CALL lim_rhg_2( i_j1, i_jpj ) 
     111               ELSE                                ! EVP rheology 
     112                  i_j1 = MAX( 1, i_j1-2 ) 
     113                  CALL lim_rhg( i_j1, i_jpj ) 
     114               ENDIF 
    122115               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, 'ij_jpj = ', i_jpj 
    123                CALL lim_rhg( i_j1, i_jpj ) 
    124                !  
    125 #endif 
     116               ! 
    126117               ! Southern hemisphere 
    127118               i_j1  =  1  
     
    130121                  i_jpj = i_jpj - 1 
    131122               END DO 
    132 #if defined key_lim2_vp 
    133                i_jpj = MIN( jpj, i_jpj+2 ) 
    134                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    135                !  
    136                CALL lim_rhg_2( i_j1, i_jpj ) 
    137                !  
    138 #else 
    139                i_jpj = MIN( jpj, i_jpj+1 ) 
     123               IF( lk_lim2_vp )   THEN             ! VP  rheology 
     124                  i_jpj = MIN( jpj, i_jpj+2 ) 
     125                  CALL lim_rhg_2( i_j1, i_jpj ) 
     126               ELSE                                ! EVP rheology 
     127                  i_jpj = MIN( jpj, i_jpj+1 ) 
     128                  CALL lim_rhg( i_j1, i_jpj ) 
     129               ENDIF 
    140130               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, 'ij_jpj = ', i_jpj 
    141                CALL lim_rhg( i_j1, i_jpj ) 
    142                !  
    143 #endif 
    144  
     131               ! 
    145132            ELSE                                 ! local domain extends over one hemisphere only 
    146133               !                                 ! Rheology is computed only over the ice cover 
     
    156143                  i_jpj = i_jpj - 1 
    157144               END DO 
    158                i_jpj = MIN( jpj, i_jpj+2) 
    159      
     145               i_jpj = MIN( jpj, i_jpj+2 ) 
     146               !  
     147               IF( lk_lim2_vp )   THEN             ! VP  rheology 
     148                  i_jpj = MIN( jpj, i_jpj+2 ) 
     149                  CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology 
     150               ELSE                                ! EVP rheology 
     151                  i_j1  = MAX( 1  , i_j1-2  ) 
     152                  i_jpj = MIN( jpj, i_jpj+1 ) 
     153                  CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology 
     154               ENDIF 
    160155               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    161                !  
    162 #if defined key_lim2_vp 
    163                i_jpj = MIN( jpj, i_jpj+2) 
    164                CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology 
    165 #else 
    166                i_j1 = MAX( 1, i_j1-2 ) 
    167                i_jpj = MIN( jpj, i_jpj+1) 
    168                CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology 
    169 #endif 
    170156               ! 
    171157            ENDIF 
     
    177163         ! computation of friction velocity 
    178164         ! -------------------------------- 
    179          SELECT CASE( cl_grid ) 
    180          CASE( 'C' )       ! C-grid ice dynamics (EVP) 
    181                            ! ice-ocean & ice velocity at  ocean velocity points 
    182             zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)   
     165         SELECT CASE( cp_ice_msh )           ! ice-ocean relative velocity at u- & v-pts 
     166         CASE( 'C' )                               ! EVP : C-grid ice dynamics 
     167            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)           ! ice-ocean & ice velocity at ocean velocity points 
    183168            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 
    184             ! 
    185          CASE( 'B' )       ! B-grid ice dynamics (VP)  
    186                            ! ice-ocean velocity at U & V-points (u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points) 
    187             DO jj = 1, jpjm1 
    188                DO ji = 1, jpim1   ! NO vector opt. 
    189                   zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
    190                   zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
     169         CASE( 'I' )                               ! VP  : B-grid ice dynamics (I-point)  
     170            DO jj = 1, jpjm1                               ! u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points 
     171               DO ji = 1, jpim1   ! NO vector opt.         ! 
     172                  zu_io(ji,jj) = 0.5_wp * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
     173                  zv_io(ji,jj) = 0.5_wp * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
    191174               END DO 
    192175            END DO 
     
    194177 
    195178         ! frictional velocity at T-point 
     179         zcoef = 0.5_wp * cw 
    196180         DO jj = 2, jpjm1 
    197181            DO ji = 2, jpim1   ! NO vector opt. because of zu_io 
    198                ust2s(ji,jj) = 0.5 * cw                                                          & 
    199                   &         * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
    200                   &            + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj) 
     182               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
     183                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj) 
    201184            END DO 
    202185         END DO 
     
    207190         DO jj = 2, jpjm1 
    208191            DO ji = fs_2, fs_jpim1   ! vector opt. 
    209                ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
    210                   &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) 
     192               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
     193                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj) 
    211194            END DO 
    212195         END DO 
     
    217200      ! 
    218201      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    219  
     202      ! 
    220203   END SUBROUTINE lim_dyn_2 
    221204 
     
    265248         WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp 
    266249      ENDIF 
     250      ! 
     251      IF( angvg /= 0._wp .AND. .NOT.lk_lim2_vp ) THEN 
     252         CALL ctl_warn( 'lim_dyn_init_2: turning angle for oceanic stress not properly coded for EVP ',   & 
     253            &           '(see limsbc_2 module). We force  angvg = 0._wp'  ) 
     254         angvg = 0._wp 
     255      ENDIF 
    267256 
    268257      !  Initialization 
Note: See TracChangeset for help on using the changeset viewer.