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 2222 for branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/limdyn_2.F90 – NEMO

Ignore:
Timestamp:
2010-10-12T15:53:09+02:00 (14 years ago)
Author:
cbricaud
Message:

commit changes from devmercator2010_1

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/limdyn_2.F90

    r2219 r2222  
    88   !!             2.0  !  03-08  (C. Ethe) add lim_dyn_init 
    99   !!             2.0  !  06-07  (G. Madec)  Surface module 
     10   !!             3.3  !  09-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    USE limrhg_2       ! ice rheology 
    25  
    26    USE lbclnk         ! 
    27    USE lib_mpp        ! 
    28    USE in_out_manager ! I/O manager 
    29    USE prtctl         ! Print control 
     19   USE dom_oce          ! ocean domain 
     20   USE sbc_oce          ! surface boundary condition variables 
     21   USE phycst           ! physical constant 
     22   USE ice_2            ! LIM2: ice variables 
     23   USE dom_ice_2        ! LIM2: ice domain 
     24   USE limistate_2      ! LIM2: ice initial state 
     25#if defined key_lim2_vp 
     26   USE limrhg_2         ! LIM2:  VP ice rheology 
     27#else 
     28   USE limrhg           ! LIM : EVP ice rheology 
     29#endif 
     30   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     31   USE lib_mpp          ! MPP library 
     32   USE in_out_manager   ! I/O manager 
     33   USE prtctl           ! Print control 
    3034 
    3135   IMPLICIT NONE 
    3236   PRIVATE 
    3337 
    34    PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim 
    35  
    36    !! * Module variables 
    37    REAL(wp)  ::  rone    = 1.e0   ! constant value 
     38   PUBLIC   lim_dyn_2   ! routine called by sbc_ice_lim module 
     39 
     40   REAL(wp) ::  rone = 1.e0   ! constant value 
    3841 
    3942#  include "vectopt_loop_substitute.h90" 
    4043   !!---------------------------------------------------------------------- 
    41    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     44   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4245   !! $Id$ 
    4346   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt) 
    4447   !!---------------------------------------------------------------------- 
    45  
    4648CONTAINS 
    4749 
     
    7173      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only) 
    7274       
    73       IF( ln_limdyn ) THEN 
    74          ! 
    75          ! Mean ice and snow thicknesses.           
    76          hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:) 
     75      IF( ln_limdyn ) THEN                     ! Rheology (ice dynamics) 
     76         !                                     ! ======== 
     77         !          
     78         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)   ! Mean ice and snow thicknesses 
    7779         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:) 
    7880         ! 
    79          !                                     ! Rheology (ice dynamics) 
    80          !                                     ! ======== 
     81         !                                     ! Define the j-limits where ice rheology is computed 
    8182          
    8283         !  Define the j-limits where ice rheology is computed 
     
    8788            i_jpj = jpj 
    8889            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    89             CALL lim_rhg_2( i_j1, i_jpj ) 
     90#if defined key_lim2_vp 
     91            CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology 
     92#else 
     93            CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology 
     94#endif 
     95         ELSE                                      !==  optimization of the computational area  ==! 
     96            DO jj = 1, jpj 
     97               zind(jj) = SUM( frld (:,jj  ) )        ! = FLOAT(jpj) if ocean everywhere on a j-line 
     98               zmsk(jj) = SUM( tmask(:,jj,1) )        ! = 0          if land  everywhere on a j-line 
     99            END DO 
    90100            ! 
    91          ELSE                                 ! optimization of the computational area 
    92             ! 
    93             DO jj = 1, jpj 
    94                zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line 
    95                zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line 
    96             END DO 
    97             ! 
    98             IF( l_jeq ) THEN                     ! local domain include both hemisphere 
    99                !                                 ! Rheology is computed in each hemisphere 
    100                !                                 ! only over the ice cover latitude strip 
    101                ! Northern hemisphere 
    102                i_j1  = njeq 
     101            IF( l_jeq ) THEN                          ! local domain include both hemisphere: rheology is computed 
     102               !                                      ! in each hemisphere only over the ice cover latitude strip 
     103               i_j1  = njeq                  ! Northern hemisphere 
    103104               i_jpj = jpj 
    104105               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
    105106                  i_j1 = i_j1 + 1 
    106107               END DO 
     108#if defined key_lim2_vp 
    107109               i_j1 = MAX( 1, i_j1-1 ) 
    108110               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    109111               !  
    110112               CALL lim_rhg_2( i_j1, i_jpj ) 
    111                !  
    112                ! Southern hemisphere 
    113                i_j1  =  1  
     113#else 
     114               i_j1 = MAX( 1, i_j1-2 ) 
     115               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     116               CALL lim_rhg( i_j1, i_jpj ) 
     117#endif 
     118               i_j1  =  1                    ! Southern hemisphere 
    114119               i_jpj = njeq 
    115120               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    116121                  i_jpj = i_jpj - 1 
    117122               END DO 
     123#if defined key_lim2_vp 
    118124               i_jpj = MIN( jpj, i_jpj+2 ) 
    119125               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    120126               !  
    121127               CALL lim_rhg_2( i_j1, i_jpj ) 
     128#else 
     129               i_jpj = MIN( jpj, i_jpj+1 ) 
     130               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     131               CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt ) 
     132#endif 
    122133               !  
    123             ELSE                                 ! local domain extends over one hemisphere only 
    124                !                                 ! Rheology is computed only over the ice cover 
    125                !                                 ! latitude strip 
     134            ELSE                                      ! local domain extends over one hemisphere only: rheology is 
     135               !                                      ! computed only over the ice cover latitude strip 
    126136               i_j1  = 1 
    127137               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
     
    129139               END DO 
    130140               i_j1 = MAX( 1, i_j1-1 ) 
    131      
    132141               i_jpj  = jpj 
    133142               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    134143                  i_jpj = i_jpj - 1 
    135144               END DO 
     145               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     146#if defined key_lim2_vp 
    136147               i_jpj = MIN( jpj, i_jpj+2) 
    137      
    138                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    139                !  
    140                CALL lim_rhg_2( i_j1, i_jpj ) 
     148               CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology 
     149#else 
     150               i_j1 = MAX( 1, i_j1-2 ) 
     151               i_jpj = MIN( jpj, i_jpj+1) 
     152               CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology 
     153#endif 
    141154               ! 
    142155            ENDIF 
    143156            ! 
    144157         ENDIF 
    145  
     158         ! 
    146159         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 
     160         !  
    147161          
    148          ! computation of friction velocity 
    149          ! -------------------------------- 
    150          ! ice-ocean velocity at U & V-points (u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points) 
    151           
    152          DO jj = 1, jpjm1 
    153             DO ji = 1, jpim1   ! NO vector opt. 
    154                zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
    155                zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
    156             END DO 
    157          END DO 
    158          ! frictional velocity at T-point 
    159          DO jj = 2, jpjm1 
     162         !                                     ! friction velocity 
     163         !                                     ! ================= 
     164         SELECT CASE( cl_grid ) 
     165         CASE( 'C' )                          ! C-grid ice dynamics (EVP) 
     166            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)      ! ice-ocean & ice velocity at ocean velocity points  
     167            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 
     168            ! 
     169         CASE( 'B' )                          ! B-grid ice dynamics (VP) 
     170            DO jj = 1, jpjm1                          ! ice velocity at I-point, ice-ocean velocity at ocean points 
     171               DO ji = 1, jpim1   ! NO vector opt. 
     172                  zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
     173                  zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
     174               END DO 
     175            END DO 
     176         END SELECT 
     177         ! 
     178         DO jj = 2, jpjm1                     ! frictional velocity at T-point 
    160179            DO ji = 2, jpim1   ! NO vector opt. because of zu_io 
    161180               ust2s(ji,jj) = 0.5 * cw                                                          & 
     
    165184         END DO 
    166185         ! 
    167       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
    168          ! 
     186      ELSE                                     ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
     187         !                                     ! =============== 
    169188         zcoef = SQRT( 0.5 ) / rau0 
    170189         DO jj = 2, jpjm1 
     
    177196      ENDIF 
    178197      ! 
    179       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point 
     198      CALL lbc_lnk( ust2s, 'T',  1. )   ! lateral boundary condition 
    180199      ! 
    181200      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    182  
     201      ! 
    183202   END SUBROUTINE lim_dyn_2 
    184203 
     
    188207      !!                  ***  ROUTINE lim_dyn_init_2  *** 
    189208      !! 
    190       !! ** Purpose :   Physical constants and parameters linked to the ice 
    191       !!              dynamics 
    192       !! 
    193       !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic 
    194       !!              parameter values 
     209      !! ** Purpose :   initialisation of  the ice dynamics variables 
     210      !!               
     211      !! ** Method  :   Read the namicedyn namelist and check their values 
    195212      !! 
    196213      !! ** input   :   Namelist namicedyn 
    197214      !!------------------------------------------------------------------- 
    198       NAMELIST/namicedyn/ epsd, alpha,     & 
    199          &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
    200          &                c_rhg, etamn, creepl, ecc, ahi0 
    201       !!------------------------------------------------------------------- 
    202  
     215      NAMELIST/namicedyn/ epsd, alpha, dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
     216         &                c_rhg, etamn, creepl, ecc, ahi0, nevp, telast, alphaevp 
     217      !!------------------------------------------------------------------- 
     218      ! 
    203219      REWIND ( numnam_ice )                       ! Read Namelist namicedyn 
    204220      READ   ( numnam_ice  , namicedyn ) 
    205  
     221      ! 
    206222      IF(lwp) THEN                                ! Control print 
    207223         WRITE(numout,*) 
     
    223239         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc 
    224240         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0 
     241         WRITE(numout,*) '       number of iterations for subcycling              nevp   = ', nevp 
     242         WRITE(numout,*) '       timescale for elastic waves                      telast = ', telast 
     243         WRITE(numout,*) '       coefficient for the solution of int. stresses  alphaevp = ', alphaevp 
    225244      ENDIF 
    226  
     245      ! 
    227246      !  Initialization 
    228247      usecc2 = 1.0 / ( ecc * ecc ) 
     
    240259#else 
    241260   !!---------------------------------------------------------------------- 
    242    !!   Default option          Empty module       NO LIM 2.0 sea-ice model 
     261   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    243262   !!---------------------------------------------------------------------- 
    244263CONTAINS 
    245    SUBROUTINE lim_dyn_2         ! Empty routine 
     264   SUBROUTINE lim_dyn_2         ! Dummy routine 
    246265   END SUBROUTINE lim_dyn_2 
    247266#endif  
Note: See TracChangeset for help on using the changeset viewer.