Changeset 1922


Ignore:
Timestamp:
2010-06-09T10:24:58+02:00 (10 years ago)
Author:
rblod
Message:

Correct alternate direction switch frenquence in LIM2 and LIM3, see ticket #625 and #605

Location:
trunk/NEMO
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/LIM_SRC_2/limtrp_2.F90

    r1715 r1922  
    44   !! LIM 2.0 transport ice model : sea-ice advection/diffusion 
    55   !!====================================================================== 
     6   !! History :  LIM  !  2000-01 (UCL)  Original code 
     7   !!            2.0  !  2001-05 (G. Madec, R. Hordoir) opa norm 
     8   !!             -   !  2004-01 (G. Madec, C. Ethe)  F90, mpp 
     9   !!---------------------------------------------------------------------- 
    610#if defined key_lim2 
    711   !!---------------------------------------------------------------------- 
     
    1115   !!   lim_trp_init_2 : initialization and namelist read 
    1216   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    14    USE phycst 
    15    USE dom_oce 
     17   USE phycst          ! physical constant 
     18   USE sbc_oce         ! ocean surface boundary condition 
     19   USE dom_oce         ! ocean domain 
    1620   USE in_out_manager  ! I/O manager 
    17    USE dom_ice_2 
    18    USE ice_2 
    19    USE limistate_2 
    20    USE limadv_2 
    21    USE limhdf_2 
    22    USE lbclnk 
    23    USE lib_mpp 
     21   USE dom_ice_2       ! LIM-2 domain 
     22   USE ice_2           ! LIM-2 variables 
     23   USE limistate_2     ! LIM-2 initial state 
     24   USE limadv_2        ! LIM-2 advection 
     25   USE limhdf_2        ! LIM-2 horizontal diffusion 
     26   USE lbclnk          ! lateral boundary conditions -- MPP exchanges 
     27   USE lib_mpp         ! MPP library 
    2428 
    2529   IMPLICIT NONE 
    2630   PRIVATE 
    2731 
    28    !! * Routine accessibility 
    29    PUBLIC lim_trp_2     ! called by sbc_ice_lim_2 
    30  
    31    !! * Shared module variables 
    32    REAL(wp), PUBLIC  ::   &  !: 
    33       bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip) 
    34  
    35    !! * Module variables 
     32   PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2 
     33 
     34   REAL(wp), PUBLIC  ::   bound  = 0.e0   !: boundary condit. (0.0 no-slip, 1.0 free-slip) 
     35 
    3636   REAL(wp)  ::           &  ! constant values 
    3737      epsi06 = 1.e-06  ,  & 
     
    4444#  include "vectopt_loop_substitute.h90" 
    4545   !!---------------------------------------------------------------------- 
    46    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
     46   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2010)  
    4747   !! $Id$ 
    48    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4949   !!---------------------------------------------------------------------- 
    5050 
     
    6262      !! 
    6363      !! ** action : 
    64       !! 
    65       !! History : 
    66       !!   1.0  !  00-01 (LIM)  Original code 
    67       !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
    68       !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp 
    6964      !!--------------------------------------------------------------------- 
    7065      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    71  
    72       INTEGER  ::   ji, jj, jk,   &  ! dummy loop indices 
    73          &          initad           ! number of sub-timestep for the advection 
    74  
    75       REAL(wp) ::  &                               
    76          zindb  ,  & 
    77          zacrith, & 
    78          zindsn , & 
    79          zindic , & 
    80          zusvosn, & 
    81          zusvoic, & 
    82          zignm  , & 
    83          zindhe , & 
    84          zvbord , & 
    85          zcfl   , & 
    86          zusnit , & 
    87          zrtt, ztsn, ztic1, ztic2 
    88  
    89       REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace 
    90          zui_u , zvi_v , zsm   ,         & 
    91          zs0ice, zs0sn , zs0a  ,         & 
    92          zs0c0 , zs0c1 , zs0c2 ,         & 
    93          zs0st 
     66      !! 
     67      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     68      INTEGER  ::   initad       ! number of sub-timestep for the advection 
     69      REAL(wp) ::   zindb  , zindsn , zindic, zacrith   ! local scalars 
     70      REAL(wp) ::   zusvosn, zusvoic, zignm , zindhe    !   -      - 
     71      REAL(wp) ::   zvbord , zcfl   , zusnit            !   -      - 
     72      REAL(wp) ::   zrtt   , ztsn   , ztic1 , ztic2     !   -      - 
     73      REAL(wp), DIMENSION(jpi,jpj)  ::   zui_u , zvi_v , zsm             ! 2D workspace 
     74      REAL(wp), DIMENSION(jpi,jpj)  ::   zs0ice, zs0sn , zs0a            !  -      - 
     75      REAL(wp), DIMENSION(jpi,jpj)  ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      - 
    9476      !--------------------------------------------------------------------- 
    9577 
     
    10587         ! ice velocities at ocean U- and V-points (zui_u,zvi_v) 
    10688         ! --------------------------------------- 
    107          ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.         
    108          zvbord = 1.0 + ( 1.0 - bound ) 
     89         zvbord = 1.0 + ( 1.0 - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions         
    10990         DO jj = 1, jpjm1 
    11091            DO ji = 1, jpim1   ! NO vector opt. 
     
    11394            END DO 
    11495         END DO 
    115          ! Lateral boundary conditions on zui_u, zvi_v 
    116          CALL lbc_lnk( zui_u, 'U', -1. ) 
    117          CALL lbc_lnk( zvi_v, 'V', -1. ) 
     96         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )         ! Lateral boundary conditions 
     97 
    11898 
    11999         ! CFL test for stability 
     
    122102         zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) ) 
    123103         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) ) 
    124  
    125          IF (lk_mpp ) CALL mpp_max(zcfl) 
    126  
    127          IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl 
     104         ! 
     105         IF(lk_mpp)   CALL mpp_max( zcfl ) 
     106         ! 
     107         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl 
    128108 
    129109         ! content of properties 
    130110         ! --------------------- 
    131111         zs0sn (:,:) =  hsnm(:,:) * area(:,:)                 ! Snow volume. 
    132          zs0ice(:,:) =  hicm (:,:) * area(:,:)                ! Ice volume. 
    133          zs0a  (:,:) =  ( 1.0 - frld(:,:) ) * area(:,:)       ! Surface covered by ice. 
    134          zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn(:,:)   ! Heat content of the snow layer. 
     112         zs0ice(:,:) =  hicm(:,:) * area(:,:)                 ! Ice volume. 
     113         zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice. 
     114         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn (:,:)  ! Heat content of the snow layer. 
    135115         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer. 
    136116         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer. 
    137          zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a(:,:)    ! Heat reservoir for brine pockets. 
     117         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)  ! Heat reservoir for brine pockets. 
    138118          
    139119  
    140          ! Advection  
     120         ! Advection (Prather scheme) 
    141121         ! --------- 
    142          ! If ice drift field is too fast, use an appropriate time step for advection.          
    143          initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) 
    144          zusnit = 1.0 / REAL( initad )  
    145           
    146          IF ( MOD( nday , 2 ) == 0) THEN 
    147             DO jk = 1,initad 
     122         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )   ! If ice drift field is too fast,           
     123         zusnit = 1.0 / REAL( initad )                              ! split the ice time step in two 
     124         ! 
     125         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==! 
     126            DO jk = 1, initad 
    148127               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) 
    149128               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) 
     
    161140               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  ) 
    162141            END DO 
    163          ELSE 
     142         ELSE                                                 !==  even ice time step:  adv_x then adv_y  ==! 
    164143            DO jk = 1, initad 
    165144               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) 
     
    182161         ! recover the properties from their contents 
    183162         ! ------------------------------------------ 
     163!!gm Define in limmsh one for all area = 1 /area  (CPU time saved !) 
    184164         zs0ice(:,:) = zs0ice(:,:) / area(:,:) 
    185165         zs0sn (:,:) = zs0sn (:,:) / area(:,:) 
     
    205185            END DO 
    206186         END DO 
     187!!gm more readable coding: (and avoid an error in F90 with sign of zero) 
     188!        DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row 
     189!           DO ji = 1 , fs_jpim1   ! vector opt. 
     190!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0.e0 
     191!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0.e0 
     192!           END DO 
     193!        END DO 
     194!!gm end 
    207195 
    208196         ! diffusion 
     
    216204         CALL lim_hdf_2( zs0st  ) 
    217205 
    218          zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) )    !!bug:  est-ce utile 
    219          zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) )    !!bug:  cf /area  juste apres 
    220          zs0a  (:,:) = MAX( rzero, zs0a  (:,:) * area(:,:) )    !! suppression des 2 change le resultat... 
    221          zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) ) 
     206!!gm see comment this can be skipped 
     207         zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) )    !!bug:  useless 
     208         zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) )    !!bug:  cf /area  just below 
     209         zs0a  (:,:) = MAX( rzero, zs0a  (:,:) * area(:,:) )    !! caution: the suppression of the 2 changes  
     210         zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) )    !! the last digit of the results 
    222211         zs0c1 (:,:) = MAX( rzero, zs0c1 (:,:) * area(:,:) ) 
    223212         zs0c2 (:,:) = MAX( rzero, zs0c2 (:,:) * area(:,:) ) 
     
    225214 
    226215 
    227          ! -------------------------------------------------------------------! 
    228          !   Up-dating and limitation of sea ice properties after transport   ! 
    229          ! -------------------------------------------------------------------! 
    230  
    231          ! Up-dating and limitation of sea ice properties after transport. 
     216         !-------------------------------------------------------------------! 
     217         !   Updating and limitation of sea ice properties after transport   ! 
     218         !-------------------------------------------------------------------! 
    232219         DO jj = 1, jpj 
    233 !!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq! 
    234220            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH 
    235221            DO ji = 1, jpi 
    236  
     222               ! 
    237223               ! Recover mean values over the grid squares. 
    238224               zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) ) 
     
    272258            END DO 
    273259         END DO 
    274           
     260         ! 
    275261      ENDIF 
    276        
     262      ! 
    277263   END SUBROUTINE lim_trp_2 
    278264 
     
    284270      !! ** Purpose :   initialization of ice advection parameters 
    285271      !! 
    286       !! ** Method  : Read the namicetrp namelist and check the parameter  
    287       !!       values called at the first timestep (nit000) 
     272      !! ** Method  :   Read the namicetrp namelist and check the parameter  
     273      !!              values called at the first timestep (nit000) 
    288274      !! 
    289275      !! ** input   :   Namelist namicetrp 
    290       !! 
    291       !! history : 
    292       !!   2.0  !  03-08 (C. Ethe)  Original code 
    293276      !!------------------------------------------------------------------- 
    294277      NAMELIST/namicetrp/ bound 
    295278      !!------------------------------------------------------------------- 
    296  
    297       ! Read Namelist namicetrp 
    298       REWIND ( numnam_ice ) 
     279      ! 
     280      REWIND ( numnam_ice )      ! Read Namelist namicetrp 
    299281      READ   ( numnam_ice  , namicetrp ) 
    300282      IF(lwp) THEN 
     
    304286         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound 
    305287      ENDIF 
    306              
     288      ! 
    307289   END SUBROUTINE lim_trp_init_2 
    308290 
  • trunk/NEMO/LIM_SRC_3/limtrp.F90

    r1715 r1922  
    204204         zusnit = 1.0 / REAL( initad )  
    205205 
    206          IF ( MOD( nday , 2 ) == 0) THEN 
     206         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==! 
    207207            DO jk = 1,initad 
    208208               !--- ice open water area 
Note: See TracChangeset for help on using the changeset viewer.