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 2135 for branches/devmercator2010_1/NEMO/LIM_SRC_2 – NEMO

Ignore:
Timestamp:
2010-09-29T19:31:33+02:00 (14 years ago)
Author:
cbricaud
Message:

add changes from branch dev_1784_EVP

Location:
branches/devmercator2010_1/NEMO/LIM_SRC_2
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/dom_ice_2.F90

    r1228 r2135  
    44   !! LIM 2.0 Sea Ice :   Domain  variables 
    55   !!====================================================================== 
    6    !! History :   2.0  !  03-08  (C. Ethe)  Free form and module 
     6   !! History :   1.0  !  2003-08  (C. Ethe)  Free form and module 
     7   !!             3.3  !  2009-05 (G.Garric, C. Bricaud) addition of lim2_evp case 
    78   !!---------------------------------------------------------------------- 
    8 #if defined key_lim2 
     9#if defined   key_lim2 
    910   !!---------------------------------------------------------------------- 
    10    !!   LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    11    !! $Id$ 
    12    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     11   !!   'key_lim2'                                       LIM2 sea-ice model 
    1312   !!---------------------------------------------------------------------- 
    1413   USE par_ice_2 
    … …  
    2019 
    2120   INTEGER, PUBLIC ::   njeq , njeqm1          !: j-index of the equator if it is inside the domain 
    22       !                                        !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
     21   !                                           !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
    2322 
    24    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor,   &  !: coriolis factor and coeficient 
    25       &                                              covrai ,         &  !: sine of geographic latitude 
    26       &                                              area   ,         &  !: surface of grid cell  
    27       &                                              tms    , tmu        !: temperature and velocity points masks 
    28    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght   ,         &  !: weight of the 4 neighbours to compute averages 
    29       &                                              akappa , bkappa     !: first and third group of metric coefficients 
    30    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd   !: second group of metric coefficients 
     23   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor     !: coriolis factor and coeficient 
     24   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   covrai            !: sine of geographic latitude 
     25   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   area              !: surface of grid cell  
     26   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tms    , tmu      !: temperature and velocity points masks 
     27   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght              !: weight of the 4 neighbours to compute averages 
    3128 
     29# if defined key_lim2_vp 
     30   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   akappa , bkappa   !: first and third group of metric coefficients 
     31   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd            !: second group of metric coefficients 
     32# else 
     33   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmv    , tmf      !: y-velocity and F-points masks 
     34   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmi               !: ice mask: =1 if ice thick > 0 
     35# endif 
     36 
     37#else 
     38   !!---------------------------------------------------------------------- 
     39   !!   Default option          Empty module         NO LIM-2 sea-ice model 
     40   !!---------------------------------------------------------------------- 
     41#endif 
     42 
     43   !!---------------------------------------------------------------------- 
     44   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
     45   !! $Id$ 
     46   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3247   !!====================================================================== 
    33 #endif 
    3448END MODULE dom_ice_2 

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/ice_2.F90

    r1756 r2135  
    44   !! Sea Ice physics:  diagnostics variables of ice defined in memory 
    55   !!===================================================================== 
    6    !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module 
     6   !! History :  2.0  !  2003-08  (C. Ethe)  F90: Free form and module 
     7   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_lim2 
    … …  
    1011   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1112   !!---------------------------------------------------------------------- 
    12    !!  LIM 2.0, UCL-LOCEAN-IPSL (2006) 
    13    !! $Id$ 
    14    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    15    !!---------------------------------------------------------------------- 
    16    !! * Modules used 
    1713   USE par_ice_2          ! LIM sea-ice parameters 
    1814 
    … …  
    2016   PRIVATE 
    2117 
    22    !!* Share parameters namelist (namicerun read in iceini) * 
     18   !                                                                     !!* namicerun namelist * 
    2319   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input) 
    2420   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output) 
    2521   LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F) 
    2622   LOGICAL               , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping 
    27    REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of temp. in snow (0) or not (9999) 
    28    REAL(wp)              , PUBLIC ::   hicdif        = 0.e0               !: computation of temp. in ice (0) or not (9999) 
     23   LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag grid points output (T) or not (F) 
     24   REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of snow temp. (0) or not (9999) 
     25   REAL(wp)              , PUBLIC ::   hicdif        = 0.e0               !: computation of ice  temp. (0) or not (9999) 
    2926   REAL(wp), DIMENSION(2), PUBLIC ::   acrit = (/ 1.e-06 , 1.e-06 /)      !: minimum fraction for leads in  
    3027   !                                                                      !: north and south hemisphere 
    31    !!* ice-dynamic namelist (namicedyn) * 
     28   !                                       !!* ice-dynamic namelist (namicedyn) * 
    3229   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation 
    3330   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation 
    … …  
    4643   REAL(wp), PUBLIC ::   ecc    = 2.e0      !: eccentricity of the elliptical yield curve 
    4744   REAL(wp), PUBLIC ::   ahi0   = 350.e0    !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    48  
     45   INTEGER , PUBLIC ::   nevp   =   360     !: number of EVP subcycling iterations 
     46   INTEGER , PUBLIC ::   telast =  3600     !: timescale for EVP elastic waves 
     47   REAL(wp), PUBLIC ::   alphaevp = 1.e0    !: coefficient for the solution of EVP int. stresses 
    4948   REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc ) 
    5049   REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw 
    … …  
    5453   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
    5554   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
    56    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses 
    57    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity 
     55   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s         !: friction velocity 
    5856 
    59    !!* diagnostic quantities 
    60    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin) 
     57# if defined key_lim2_vp 
     58   !                                                        !!* VP rheology * 
     59   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm     !: mean snow and ice thicknesses 
     60   CHARACTER(len=1), PUBLIC :: cl_grid =   'B'                   !: type of grid used in ice dynamics, 'C' or 'B' 
     61# else 
     62   !                                                        !!* EVP rheology * 
     63   CHARACTER(len=1), PUBLIC             ::   cl_grid = 'C'   !: type of grid used in ice dynamics, 'C' or 'B' 
     64   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress1_i       !: first stress tensor element        
     65   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress2_i       !: second stress tensor element 
     66   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress12_i      !: diagonal stress tensor element 
     67   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   delta_i         !: rheology delta factor (see Flato and Hibler 95) [s-1] 
     68   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   divu_i          !: Divergence of the velocity field [s-1] -> limrhg.F90 
     69   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   shear_i         !: Shear of the velocity field [s-1] -> limrhg.F90 
     70   ! 
     71   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)          :: at_i          !: 
     72   REAL(wp), PUBLIC, DIMENSION(:,:)    , POINTER :: vt_s ,vt_i    !: mean snow and ice thicknesses 
     73   REAL(wp), PUBLIC, DIMENSION(jpi,jpj), TARGET  :: hsnm , hicm   !: target vt_s ,vt_i pointers  
     74#endif 
     75 
     76   !                                                      !!* diagnostic quantities 
     77   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvosif       !: ice volume change at ice surface (only used for outputs) 
     78   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvobif       !: ice volume change at ice bottom  (only used for outputs) 
     79   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdvolif       !: Total   ice volume change (only used for outputs) 
     80   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvonif       !: Lateral ice volume change (only used for outputs) 
     81   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature [Kelvin] 
    6182   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS 
    6283   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness 
    … …  
    7192   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass 
    7293   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean) 
    73    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing  
     94   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring  ocean surface layer to freezing  
    7495   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean 
    7596   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2) 
    7697   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget 
    7798   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice 
    78    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max heat contained in brine pockets (?) 
     99   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max brine pockets heat content 
    79100   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?) 
    80101   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?) 
    81    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy in case of toral lateral ablation (?) 
     102   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy (total lateral ablation case) 
    82103   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice 
    83104 
    84    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice   !: two components of the ice   velocity at I-point (m/s) 
    85    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce   !: two components of the ocean velocity at I-point (m/s) 
     105   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice   !: two components of the ice   velocity at I-point [m/s] 
     106   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce   !: two components of the ocean velocity at I-point [m/s] 
    86107 
    87108   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer 
    88109 
    89    !!* moment used in the advection scheme 
     110   !                                                                               !!* moment used in the advection 
    90111   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume 
    91112   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume 
    … …  
    98119#else 
    99120   !!---------------------------------------------------------------------- 
    100    !!   Default option         Empty module        NO LIM 2.0 sea-ice model 
     121   !!   Default option          Empty module         NO LIM-2 sea-ice model 
    101122   !!---------------------------------------------------------------------- 
    102123#endif 
    103124 
     125   !!---------------------------------------------------------------------- 
     126   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
     127   !! $Id$ 
     128   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    104129   !!====================================================================== 
    105130END MODULE ice_2 

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/iceini_2.F90

    r1581 r2135  
    66   !! History :   1.0  !  02-08  (G. Madec)  F90: Free form and modules 
    77   !!             2.0  !  03-08  (C. Ethe)  add ice_run 
     8   !!             3.3  !  09-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_lim2 
    … …  
    1112   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1213   !!---------------------------------------------------------------------- 
    13    !!---------------------------------------------------------------------- 
    1414   !!   ice_init_2       : sea-ice model initialization 
    1515   !!   ice_run_2        : Definition some run parameter for ice model 
    1616   !!---------------------------------------------------------------------- 
    17    USE dom_oce 
    18    USE dom_ice_2 
    19    USE sbc_oce         ! surface boundary condition: ocean 
    20    USE sbc_ice         ! surface boundary condition: ice 
    21    USE phycst          ! Define parameters for the routines 
    22    USE ice_2 
    23    USE limmsh_2 
    24    USE limistate_2 
    25    USE limrst_2    
    26    USE in_out_manager 
     17   USE dom_oce          ! ocean domain 
     18   USE dom_ice_2        ! LIM2: ice domain 
     19   USE sbc_oce          ! surface boundary condition: ocean 
     20   USE sbc_ice          ! surface boundary condition: ice 
     21   USE phycst           ! Define parameters for the routines 
     22   USE ice_2            ! LIM2: ice variable 
     23   USE limmsh_2         ! LIM2: mesh 
     24   USE limistate_2      ! LIM2: initial state 
     25   USE limrst_2         ! LIM2: restart 
     26   USE in_out_manager   ! I/O manager 
    2727       
    2828   IMPLICIT NONE 
    2929   PRIVATE 
    3030 
    31    PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90 
     31   PUBLIC   ice_init_2   ! called by sbcice_lim_2.F90 
     32 
     33   INTEGER, PUBLIC ::   numit   !: iteration number 
    3234 
    3335   !!---------------------------------------------------------------------- 
    34    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    35    !! $Id$  
    36    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     36   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
     37   !! $Id$ 
     38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3739   !!---------------------------------------------------------------------- 
    38  
    3940CONTAINS 
    4041 
    … …  
    4344      !!                  ***  ROUTINE ice_init_2  *** 
    4445      !! 
    45       !! ** purpose :    
     46      !! ** purpose :   initialisation of LIM-2 domain and variables  
    4647      !!---------------------------------------------------------------------- 
    4748      ! 
    48       ! Open the namelist file  
     49      !                               ! Open the namelist file  
    4950      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )       
    50       CALL ice_run_2                    !  read in namelist some run parameters 
     51      CALL ice_run_2                  !  read in namelist some run parameters 
    5152                  
    52       ! Louvain la Neuve Ice model 
    53       rdt_ice = nn_fsbc * rdttra(1) 
     53      rdt_ice = nn_fsbc * rdttra(1)   ! sea-ice time step 
     54      numit   = nit000 - 1               
    5455 
    5556      CALL lim_msh_2                  ! ice mesh initialization 
    5657      
    57       ! Initial sea-ice state 
    58       IF( .NOT.ln_rstart ) THEN 
    59          CALL lim_istate_2            ! start from rest: sea-ice deduced from sst 
    60       ELSE 
    61          CALL lim_rst_read_2          ! start from a restart file 
     58      !                               ! ice initialisation (start from rest or from a restart) 
     59      IF( .NOT.ln_rstart ) THEN   ;   CALL lim_istate_2      
     60      ELSE                        ;   CALL lim_rst_read_2 
    6261      ENDIF 
    6362       
    64       tn_ice(:,:,1) = sist(:,:)         ! initialisation of ice temperature    
     63      tn_ice(:,:,1) = sist(:,:)       ! initialisation of ice temperature    
    6564      fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction     
    6665      ! 
    … …  
    7574      !! 
    7675      !! ** Method  :   Read the namicerun namelist and check the parameter  
    77       !!       values called at the first timestep (nit000) 
     76      !!              values called at the first timestep (nit000) 
    7877      !! 
    7978      !! ** input   :   Namelist namicerun 
    … …  
    8281      !!------------------------------------------------------------------- 
    8382      !                     
    84       REWIND ( numnam_ice )                       ! Read Namelist namicerun  
    85       READ   ( numnam_ice , namicerun ) 
    86  
    87       IF(lwp) THEN 
     83      REWIND( numnam_ice )                   ! Read namicerun namelist 
     84      READ  ( numnam_ice , namicerun ) 
     85      ! 
     86      IF(lwp) THEN                           ! control print 
    8887         WRITE(numout,*) 
    8988         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/limdyn_2.F90

    r1694 r2135  
    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 (modipsl/doc/NEMO_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          
    82          !  Define the j-limits where ice rheology is computed 
    83          ! --------------------------------------------------- 
    84           
    85          IF( lk_mpp .OR. nbit_cmp == 1 ) THEN                    ! mpp: compute over the whole domain 
     83         IF( lk_mpp .OR. nbit_cmp == 1 ) THEN      !==  mpp: compute over the whole domain  ==! 
    8684            i_j1 = 1    
    8785            i_jpj = jpj 
    8886            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 ) 
     87#if defined key_lim2_vp 
     88            CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology 
     89#else 
     90            CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology 
     91#endif 
     92         ELSE                                      !==  optimization of the computational area  ==! 
     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 
    9097            ! 
    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 
     98            IF( l_jeq ) THEN                          ! local domain include both hemisphere: rheology is computed 
     99               !                                      ! in each hemisphere only over the ice cover latitude strip 
     100               i_j1  = njeq                  ! Northern hemisphere 
    103101               i_jpj = jpj 
    104102               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
    105103                  i_j1 = i_j1 + 1 
    106104               END DO 
     105#if defined key_lim2_vp 
    107106               i_j1 = MAX( 1, i_j1-1 ) 
    108107               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    109108               !  
    110109               CALL lim_rhg_2( i_j1, i_jpj ) 
    111                !  
    112                ! Southern hemisphere 
    113                i_j1  =  1  
     110#else 
     111               i_j1 = MAX( 1, i_j1-2 ) 
     112               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     113               CALL lim_rhg( i_j1, i_jpj ) 
     114#endif 
     115               i_j1  =  1                    ! Southern hemisphere 
    114116               i_jpj = njeq 
    115117               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    116118                  i_jpj = i_jpj - 1 
    117119               END DO 
     120#if defined key_lim2_vp 
    118121               i_jpj = MIN( jpj, i_jpj+2 ) 
    119122               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    120123               !  
    121124               CALL lim_rhg_2( i_j1, i_jpj ) 
     125#else 
     126               i_jpj = MIN( jpj, i_jpj+1 ) 
     127               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     128               CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt ) 
     129#endif 
    122130               !  
    123             ELSE                                 ! local domain extends over one hemisphere only 
    124                !                                 ! Rheology is computed only over the ice cover 
    125                !                                 ! latitude strip 
     131            ELSE                                      ! local domain extends over one hemisphere only: rheology is 
     132               !                                      ! computed only over the ice cover latitude strip 
    126133               i_j1  = 1 
    127134               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
    … …  
    129136               END DO 
    130137               i_j1 = MAX( 1, i_j1-1 ) 
    131      
    132138               i_jpj  = jpj 
    133139               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    134140                  i_jpj = i_jpj - 1 
    135141               END DO 
     142               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     143#if defined key_lim2_vp 
    136144               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 ) 
     145               CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology 
     146#else 
     147               i_j1 = MAX( 1, i_j1-2 ) 
     148               i_jpj = MIN( jpj, i_jpj+1) 
     149               CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology 
     150#endif 
    141151               ! 
    142152            ENDIF 
    143153            ! 
    144154         ENDIF 
    145  
     155         ! 
    146156         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 
     157         !  
    147158          
    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 
     159         !                                     ! friction velocity 
     160         !                                     ! ================= 
     161         SELECT CASE( cl_grid ) 
     162         CASE( 'C' )                          ! C-grid ice dynamics (EVP) 
     163            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)      ! ice-ocean & ice velocity at ocean velocity points  
     164            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 
     165            ! 
     166         CASE( 'B' )                          ! B-grid ice dynamics (VP) 
     167            DO jj = 1, jpjm1                          ! ice velocity at I-point, ice-ocean velocity at ocean points 
     168               DO ji = 1, jpim1   ! NO vector opt. 
     169                  zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
     170                  zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
     171               END DO 
     172            END DO 
     173         END SELECT 
     174         ! 
     175         DO jj = 2, jpjm1                     ! frictional velocity at T-point 
    160176            DO ji = 2, jpim1   ! NO vector opt. because of zu_io 
    161177               ust2s(ji,jj) = 0.5 * cw                                                          & 
    … …  
    165181         END DO 
    166182         ! 
    167       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
    168          ! 
     183      ELSE                                     ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
     184         !                                     ! =============== 
    169185         zcoef = SQRT( 0.5 ) / rau0 
    170186         DO jj = 2, jpjm1 
    … …  
    177193      ENDIF 
    178194      ! 
    179       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point 
     195      CALL lbc_lnk( ust2s, 'T',  1. )   ! lateral boundary condition 
    180196      ! 
    181197      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    182  
     198      ! 
    183199   END SUBROUTINE lim_dyn_2 
    184200 
    … …  
    188204      !!                  ***  ROUTINE lim_dyn_init_2  *** 
    189205      !! 
    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 
     206      !! ** Purpose :   initialisation of  the ice dynamics variables 
     207      !!               
     208      !! ** Method  :   Read the namicedyn namelist and check their values 
    195209      !! 
    196210      !! ** input   :   Namelist namicedyn 
    197211      !!------------------------------------------------------------------- 
    198       NAMELIST/namicedyn/ epsd, alpha,     & 
    199          &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
    200          &                c_rhg, etamn, creepl, ecc, ahi0 
    201       !!------------------------------------------------------------------- 
    202  
     212      NAMELIST/namicedyn/ epsd, alpha, dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
     213         &                c_rhg, etamn, creepl, ecc, ahi0, nevp, telast, alphaevp 
     214      !!------------------------------------------------------------------- 
     215      ! 
    203216      REWIND ( numnam_ice )                       ! Read Namelist namicedyn 
    204217      READ   ( numnam_ice  , namicedyn ) 
    205  
     218      ! 
    206219      IF(lwp) THEN                                ! Control print 
    207220         WRITE(numout,*) 
    … …  
    223236         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc 
    224237         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0 
     238         WRITE(numout,*) '       number of iterations for subcycling              nevp   = ', nevp 
     239         WRITE(numout,*) '       timescale for elastic waves                      telast = ', telast 
     240         WRITE(numout,*) '       coefficient for the solution of int. stresses  alphaevp = ', alphaevp 
    225241      ENDIF 
    226  
     242      ! 
    227243      !  Initialization 
    228244      usecc2 = 1.0 / ( ecc * ecc ) 
    … …  
    240256#else 
    241257   !!---------------------------------------------------------------------- 
    242    !!   Default option          Empty module       NO LIM 2.0 sea-ice model 
     258   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    243259   !!---------------------------------------------------------------------- 
    244260CONTAINS 
    245    SUBROUTINE lim_dyn_2         ! Empty routine 
     261   SUBROUTINE lim_dyn_2         ! Dummy routine 
    246262   END SUBROUTINE lim_dyn_2 
    247263#endif  

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/limmsh_2.F90

    r1694 r2135  
    44   !! LIM 2.0 ice model :   definition of the ice mesh parameters 
    55   !!====================================================================== 
    6 #if defined key_lim2 
     6   !! History :  LIM  ! 2001-04  (Louvain-la-Neuve) Original code 
     7   !!            1.0  ! 2002-08  (C. Ethe, G. Madec) 
     8   !!            3.3  ! 2009-05  (G. Garric, C. Bricaud) addition of the lim2_evp case 
     9   !!---------------------------------------------------------------------- 
     10#if defined   key_lim2 
    711   !!---------------------------------------------------------------------- 
    812   !!   'key_lim2'                                     LIM 2.0sea-ice model 
    … …  
    1014   !!   lim_msh_2   : definition of the ice mesh 
    1115   !!---------------------------------------------------------------------- 
    12    !! * Modules used 
    13    USE phycst 
    14    USE dom_oce 
    15    USE dom_ice_2 
    16    USE lbclnk 
    17    USE in_out_manager 
     16   USE phycst           ! physical constants 
     17   USE dom_oce          ! ocean domain 
     18   USE dom_ice_2        ! LIM2: ice domain 
     19   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     20   USE in_out_manager   ! I/O manager 
    1821 
    1922   IMPLICIT NONE 
    2023   PRIVATE 
    2124 
    22    !! * Accessibility 
    23    PUBLIC lim_msh_2      ! routine called by ice_ini_2.F90 
    24  
    25    !!---------------------------------------------------------------------- 
    26    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
     25   PUBLIC   lim_msh_2   ! routine called by ice_ini_2.F90 
     26 
     27   !!---------------------------------------------------------------------- 
     28   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    2729   !! $Id$ 
    28    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    29    !!---------------------------------------------------------------------- 
    30  
     30   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     31   !!---------------------------------------------------------------------- 
    3132CONTAINS 
    3233 
    … …  
    4243      !!              - Initialization of the ice masks (tmsk, umsk) 
    4344      !!  
    44       !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10  
    45       !! 
    46       !! ** History : 
    47       !!         original    : 01-04 (LIM) 
    48       !!         addition    : 02-08 (C. Ethe, G. Madec) 
     45      !! References : Deleersnijder et al. Ocean Modelling 100, 7-10  
    4946      !!---------------------------------------------------------------------  
    50       !! * Local variables 
    5147      INTEGER :: ji, jj      ! dummy loop indices 
    52  
    53       REAL(wp), DIMENSION(jpi,jpj) ::  & 
    54          zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction 
    55          !                   ! (resp. y direction) (defined at the center) 
    56       REAL(wp) ::         & 
    57          zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid 
    58          zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid 
    59          zusden, zusden2     ! temporary scalars 
     48      REAL(wp) ::   zusden   ! local scalars 
     49#if defined key_lim2_vp 
     50      REAL(wp) ::   zusden2           ! local scalars 
     51      REAL(wp) ::   zh1p  , zh2p      !   -      - 
     52      REAL(wp) ::   zd2d1p, zd1d2p    !   -      - 
     53      REAL(wp), DIMENSION(jpi,jpj) ::   zd2d1 , zd1d2   ! 2D workspace 
     54#endif 
    6055      !!--------------------------------------------------------------------- 
    61  
     56      ! 
    6257      IF(lwp) THEN 
    6358         WRITE(numout,*) 
    … …  
    7368      njeqm1 = njeq - 1  
    7469 
    75       fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor 
     70      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor at T-point 
    7671  
    77 !i    DO jj = 1, jpj 
    78 !i       zmsk(jj) = SUM( tmask(:,jj,:) )   ! = 0          if land  everywhere on a j-line 
    79 !!ii     write(numout,*) jj, zind(jj) 
    80 !i    END DO 
    81  
    8272      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere 
    8373         l_jeq = .TRUE. 
    … …  
    112102      !------------------- 
    113103!!ibug ??? 
    114       akappa(:,:,:,:) = 0.e0 
    115104      wght(:,:,:,:) = 0.e0 
     105      tmu(:,:)      = 0.e0 
     106#if defined key_lim2_vp  
     107      akappa(:,:,:,:)     = 0.e0 
    116108      alambd(:,:,:,:,:,:) = 0.e0 
    117       tmu(:,:) = 0.e0 
     109#else 
     110      tmv(:,:) = 0.e0 
     111      tmf(:,:) = 0.e0 
     112#endif 
    118113!!i 
    119114       
    120        
    121       ! metric coefficients for sea ice dynamic 
     115#if defined key_lim2_vp 
     116      ! metric coefficients for sea ice dynamic (VP rheology) 
    122117      !---------------------------------------- 
    123118      !                                                       ! akappa 
    … …  
    125120         zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1) 
    126121      END DO 
    127       CALL lbc_lnk( zd1d2, 'T', -1. ) 
    128  
    129122      DO ji = 2, jpi 
    130123         zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:) 
    131124      END DO 
    132       CALL lbc_lnk( zd2d1, 'T', -1. ) 
    133  
     125      CALL lbc_lnk( zd1d2, 'T', -1. )   ;   CALL lbc_lnk( zd2d1, 'T', -1. )      ! lateral boundary condition 
     126      ! 
    134127      akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) ) 
    135128      akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) ) 
    … …  
    137130      akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) ) 
    138131       
    139       !                                                      ! weights (wght) 
    140       DO jj = 2, jpj 
     132      !              
     133      DO jj = 2, jpj                                         ! weights (wght) at I-points 
    141134         DO ji = 2, jpi 
    142135            zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   & 
    … …  
    152145      CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used 
    153146      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. ) 
     147#else 
     148      ! metric coefficients for sea ice dynamic (EVP rheology) 
     149      !---------------------------------------- 
     150      DO jj = 1, jpjm1                                       ! weights (wght) at F-points 
     151         DO ji = 1, jpim1 
     152            zusden = 1. / (  ( e1t(ji+1,jj  ) + e1t(ji,jj) )   & 
     153               &           * ( e2t(ji  ,jj+1) + e2t(ji,jj) ) )  
     154            wght(ji,jj,1,1) = zusden * e1t(ji+1,jj) * e2t(ji,jj+1) 
     155            wght(ji,jj,1,2) = zusden * e1t(ji+1,jj) * e2t(ji,jj  ) 
     156            wght(ji,jj,2,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj+1) 
     157            wght(ji,jj,2,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj  ) 
     158         END DO 
     159      END DO 
     160      CALL lbc_lnk( wght(:,:,1,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,1,2), 'F', 1. )       ! lateral boundary cond.    
     161      CALL lbc_lnk( wght(:,:,2,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,2,2), 'F', 1. ) 
     162#endif 
    154163     
    155164      ! Coefficients for divergence of the stress tensor 
    156165      !------------------------------------------------- 
    157  
    158       DO jj = 2, jpj 
     166#if defined key_lim2_vp 
     167      DO jj = 2, jpj                                     ! VP rheology 
    159168         DO ji = 2, jpi   ! NO vector opt. 
    160             zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2)   & 
    161                &   + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    162                &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1)   & 
    163                &   + e1t(ji-1,jj-1) * wght(ji,jj,1,1) 
    164  
    165             zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2)   & 
    166                &   + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    167                &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1)   & 
    168                &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1) 
    169  
    170 ! better written but change the last digit and thus solver in less than 100 timestep 
    171 !           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   & 
    172 !              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)  
    173  
    174 !           zh2p  = e2t(ji-1,jj  ) * wght(ji,jj,1,2) + e2t(ji,jj  ) * wght(ji,jj,2,2)   & 
    175 !              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1) 
    176  
    177 !!ibug =0   zusden = 1.0 / ( zh1p * zh2p * 4.e0 ) 
    178             zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 ) 
     169            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2) + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     170               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1) + e1t(ji-1,jj-1) * wght(ji,jj,1,1) 
     171            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2) + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     172               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1) + e2t(ji-1,jj-1) * wght(ji,jj,1,1) 
     173            ! 
     174            zusden  = 1.e0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 ) 
    179175            zusden2 = zusden * 2.0  
    180  
    181             zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   ) 
    182             zd2d1p = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   ) 
    183  
     176            zd1d2p  = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   ) 
     177            zd2d1p  = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   ) 
     178            ! 
    184179            alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1) 
    185180            alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  ) 
    186181            alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1) 
    187182            alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  ) 
    188  
     183            ! 
    189184            alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1) 
    190185            alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  ) 
    191186            alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1) 
    192187            alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  ) 
    193  
     188            ! 
    194189            alambd(ji,jj,1,2,2,1) = zd1d2p 
    195190            alambd(ji,jj,1,2,2,2) = zd1d2p 
    196191            alambd(ji,jj,1,2,1,1) = zd1d2p 
    197192            alambd(ji,jj,1,2,1,2) = zd1d2p 
    198  
     193            ! 
    199194            alambd(ji,jj,2,1,2,1) = zd2d1p 
    200195            alambd(ji,jj,2,1,2,2) = zd2d1p 
    … …  
    203198         END DO 
    204199      END DO 
    205  
    206       CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V point 
    207       CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )      ! the value of wght at jpj is wrong 
    208       CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )      ! but it is never used 
    209       CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )      !  
    210  
    211       CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )      ! CAUTION: idem 
    212       CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )      !  
    213       CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )      ! 
    214       CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )      ! 
    215  
    216       CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )      ! CAUTION: idem 
    217       CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )      ! 
    218       CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )      ! 
    219       CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )      ! 
    220  
    221       CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )      ! CAUTION: idem 
    222       CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )      ! 
    223       CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      ! 
    224       CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      ! 
     200      ! 
     201      ! lateral boundary conditions 
     202      ! CAUTION: even with the lbc_lnk at ice U-V point, the value of wght at jpj is wrong but it is never used 
     203      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )  
     204      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )  
     205      ! 
     206      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. ) 
     207      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. ) 
     208      ! 
     209      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. ) 
     210      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. ) 
     211      ! 
     212      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. ) 
     213      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. ) 
     214#endif 
    225215             
    226216 
    227217      ! Initialization of ice masks 
    228218      !---------------------------- 
    229        
     219      !       
    230220      tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask 
    231  
    232 !i here we can use umask with a i and j shift of -1,-1 
    233       tmu(:,1) = 0.e0 
     221      ! 
     222#if defined key_lim2_vp 
     223      ! VP rheology : ice velocity point is I-point 
     224      tmu(:,1) = 0.e0              !  
    234225      tmu(1,:) = 0.e0 
    235       DO jj = 2, jpj               ! ice U.V-point: computed from ice T-point mask 
     226      DO jj = 2, jpj 
    236227         DO ji = 2, jpim1   ! NO vector opt. 
    237228            tmu(ji,jj) =  tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)             
    238229         END DO 
    239230      END DO 
    240        
    241       !--lateral boundary conditions     
    242       CALL lbc_lnk( tmu(:,:), 'I', 1. ) 
    243        
    244       ! unmasked and masked area of T-grid cell 
    245       area(:,:) = e1t(:,:) * e2t(:,:) 
    246        
     231      CALL lbc_lnk( tmu(:,:), 'I', 1. )      ! lateral boundary conditions     
     232#else 
     233      ! EVP rheology : ice velocity point are U- & V-points ; ice vorticity point is F-point 
     234      tmu(:,:) = umask(:,:,1) 
     235      tmv(:,:) = vmask(:,:,1) 
     236      tmf(:,:) = 0.e0                        ! used of fmask except its special value along the coast (rn_shlat) 
     237      WHERE( fmask(:,:,1) == 1.e0 )   tmf(:,:) = 1.e0 
     238#endif 
     239      ! 
     240      area(:,:) = e1t(:,:) * e2t(:,:)        ! unmasked and masked area of T-grid cell 
     241      ! 
    247242   END SUBROUTINE lim_msh_2 
    248243 

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/limrhg_2.F90

    r1774 r2135  
    44   !!   Ice rheology :  performs sea ice rheology 
    55   !!====================================================================== 
    6    !! History :  0.0  !  93-12  (M.A. Morales Maqueda.)  Original code 
    7    !!            1.0  !  94-12  (H. Goosse)  
    8    !!            2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp 
    9    !!            " "  !  06-08  (G. Madec)  surface module, ice-stress at I-point 
    10    !!            " "  !  09-09  (G. Madec)  Huge verctor optimisation 
    11    !!---------------------------------------------------------------------- 
    12 #if defined key_lim2 
    13    !!---------------------------------------------------------------------- 
    14    !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    15    !!---------------------------------------------------------------------- 
     6   !! History :  LIM  !  1993-12  (M.A. Morales Maqueda.)  Original code 
     7   !!            1.0  !  1994-12  (H. Goosse)  
     8   !!            2.0  !  2003-12  (C. Ethe, G. Madec)  F90, mpp 
     9   !!             -   !  2006-08  (G. Madec)  surface module, ice-stress at I-point 
     10   !!             -   !  2009-09  (G. Madec)  Huge verctor optimisation 
     11   !!            3.3  !  2009-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
     12   !!---------------------------------------------------------------------- 
     13#if defined   key_lim2   &&   defined key_lim2_vp 
     14   !!---------------------------------------------------------------------- 
     15   !!   'key_lim2'                and                 LIM 2.0 sea-ice model 
     16   !!   'key_lim2_vp'                                       VP ice rheology 
    1617   !!---------------------------------------------------------------------- 
    1718   !!   lim_rhg_2   : computes ice velocities 
    … …  
    2122   USE sbc_oce        ! surface boundary condition: ocean variables 
    2223   USE sbc_ice        ! surface boundary condition: ice variables 
    23    USE dom_ice_2      ! domaine: ice variables 
    24    USE phycst         ! physical constant 
    25    USE ice_2          ! ice variables 
    26    USE lbclnk         ! lateral boundary condition 
     24   USE dom_ice_2      ! LIM2: ice domain 
     25   USE phycst         ! physical constants 
     26   USE ice_2          ! LIM2: ice variables 
     27   USE lbclnk         ! lateral boundary condition - MPP exchanges 
    2728   USE lib_mpp        ! MPP library 
    2829   USE in_out_manager ! I/O manager 
    … …  
    3233   PRIVATE 
    3334 
    34    PUBLIC   lim_rhg_2 ! routine called by lim_dyn 
     35   PUBLIC   lim_rhg_2   ! routine called by lim_dyn 
    3536 
    3637   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero 
    … …  
    4041#  include "vectopt_loop_substitute.h90" 
    4142   !!---------------------------------------------------------------------- 
    42    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     43   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4344   !! $Id$ 
    4445   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    8990      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zzfrld, zztms 
    9091      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zi1, zi2, zmasst, zpresh 
    91  
    9292      !!------------------------------------------------------------------- 
    93  
    94 !!bug 
    95 !!    u_oce(:,:) = 0.e0 
    96 !!    v_oce(:,:) = 0.e0 
    97 !!    write(*,*) 'rhg min, max u & v', maxval(u_oce), minval(u_oce), maxval(v_oce), minval(v_oce) 
    98 !!bug 
    9993       
    10094      !  Store initial velocities 
    … …  
    106100      zztms(:,1:jpj) = tms(:,:)   ;    zzfrld(:,1:jpj) = frld(:,:) 
    107101      zu0(:,1:jpj) = u_ice(:,:)   ;    zv0(:,1:jpj) = v_ice(:,:) 
    108  
     102      ! 
    109103      zu_a(:,:)    = zu0(:,:)     ;   zv_a(:,:) = zv0(:,:) 
    110104      zu_n(:,:)    = zu0(:,:)     ;   zv_n(:,:) = zv0(:,:) 
    … …  
    126120      zviszeta(:,:) = 0.e0 
    127121      zviseta (:,:) = 0.e0 
    128  
    129 !i    zviszeta(:,0    ) = 0.e0    ;    zviseta(:,0    ) = 0.e0 
    130 !i    zviszeta(:,jpj  ) = 0.e0    ;    zviseta(:,jpj  ) = 0.e0 
    131 !i    zviszeta(:,jpj+1) = 0.e0    ;    zviseta(:,jpj+1) = 0.e0 
    132122 
    133123 
    … …  
    372362            END DO 
    373363         END DO 
    374  
    375          ! GAUSS-SEIDEL method 
     364         !  
    376365         !                                                      ! ================ ! 
    377 iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! 
     366iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! GAUSS-SEIDEL method 
    378367            !                                                   ! ================ ! 
    379368!CDIR NOVERRCHK 

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/limrst_2.F90

    r1715 r2135  
    66   !! History :  2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
    77   !!                 !  06-07  (S. Masson)  use IOM for restart read/write 
     8   !!            3.3  !  09-05  (G.Garric) addition of the lim2_evp case 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_lim2 
    … …  
    108109      CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) )  
    109110       
    110       CALL iom_rstput( iter, nitrst, numriw, 'hicif' , hicif (:,:)   )      ! prognostic variables  
    111       CALL iom_rstput( iter, nitrst, numriw, 'hsnif' , hsnif (:,:)   ) 
    112       CALL iom_rstput( iter, nitrst, numriw, 'frld'  , frld  (:,:)   ) 
    113       CALL iom_rstput( iter, nitrst, numriw, 'sist'  , sist  (:,:)   ) 
    114       CALL iom_rstput( iter, nitrst, numriw, 'tbif1' , tbif  (:,:,1) ) 
    115       CALL iom_rstput( iter, nitrst, numriw, 'tbif2' , tbif  (:,:,2) ) 
    116       CALL iom_rstput( iter, nitrst, numriw, 'tbif3' , tbif  (:,:,3) ) 
    117       CALL iom_rstput( iter, nitrst, numriw, 'u_ice' , u_ice (:,:)   ) 
    118       CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice (:,:)   ) 
    119       CALL iom_rstput( iter, nitrst, numriw, 'qstoif', qstoif(:,:)   ) 
    120       CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq (:,:)   ) 
    121       CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice (:,:)   ) 
    122       CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice (:,:)   ) 
    123       CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice(:,:)   ) 
    124       CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice(:,:)   ) 
    125       CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice(:,:)   ) 
    126       CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn  (:,:)   ) 
    127       CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn  (:,:)   ) 
    128       CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn (:,:)   ) 
    129       CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn (:,:)   ) 
    130       CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn (:,:)   ) 
    131       CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa   (:,:)   ) 
    132       CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya   (:,:)   ) 
    133       CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa  (:,:)   ) 
    134       CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya  (:,:)   ) 
    135       CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya  (:,:)   ) 
    136       CALL iom_rstput( iter, nitrst, numriw, 'sxc0'  , sxc0  (:,:)   ) 
    137       CALL iom_rstput( iter, nitrst, numriw, 'syc0'  , syc0  (:,:)   ) 
    138       CALL iom_rstput( iter, nitrst, numriw, 'sxxc0' , sxxc0 (:,:)   ) 
    139       CALL iom_rstput( iter, nitrst, numriw, 'syyc0' , syyc0 (:,:)   ) 
    140       CALL iom_rstput( iter, nitrst, numriw, 'sxyc0' , sxyc0 (:,:)   ) 
    141       CALL iom_rstput( iter, nitrst, numriw, 'sxc1'  , sxc1  (:,:)   ) 
    142       CALL iom_rstput( iter, nitrst, numriw, 'syc1'  , syc1  (:,:)   ) 
    143       CALL iom_rstput( iter, nitrst, numriw, 'sxxc1' , sxxc1 (:,:)   ) 
    144       CALL iom_rstput( iter, nitrst, numriw, 'syyc1' , syyc1 (:,:)   ) 
    145       CALL iom_rstput( iter, nitrst, numriw, 'sxyc1' , sxyc1 (:,:)   ) 
    146       CALL iom_rstput( iter, nitrst, numriw, 'sxc2'  , sxc2  (:,:)   ) 
    147       CALL iom_rstput( iter, nitrst, numriw, 'syc2'  , syc2  (:,:)   ) 
    148       CALL iom_rstput( iter, nitrst, numriw, 'sxxc2' , sxxc2 (:,:)   ) 
    149       CALL iom_rstput( iter, nitrst, numriw, 'syyc2' , syyc2 (:,:)   ) 
    150       CALL iom_rstput( iter, nitrst, numriw, 'sxyc2' , sxyc2 (:,:)   ) 
    151       CALL iom_rstput( iter, nitrst, numriw, 'sxst'  , sxst  (:,:)   ) 
    152       CALL iom_rstput( iter, nitrst, numriw, 'syst'  , syst  (:,:)   ) 
    153       CALL iom_rstput( iter, nitrst, numriw, 'sxxst' , sxxst (:,:)   ) 
    154       CALL iom_rstput( iter, nitrst, numriw, 'syyst' , syyst (:,:)   ) 
    155       CALL iom_rstput( iter, nitrst, numriw, 'sxyst' , sxyst (:,:)   ) 
     111      CALL iom_rstput( iter, nitrst, numriw, 'hicif'      , hicif     (:,:)   )      ! prognostic variables  
     112      CALL iom_rstput( iter, nitrst, numriw, 'hsnif'      , hsnif     (:,:)   ) 
     113      CALL iom_rstput( iter, nitrst, numriw, 'frld'       , frld      (:,:)   ) 
     114      CALL iom_rstput( iter, nitrst, numriw, 'sist'       , sist      (:,:)   ) 
     115      CALL iom_rstput( iter, nitrst, numriw, 'tbif1'      , tbif      (:,:,1) ) 
     116      CALL iom_rstput( iter, nitrst, numriw, 'tbif2'      , tbif      (:,:,2) ) 
     117      CALL iom_rstput( iter, nitrst, numriw, 'tbif3'      , tbif      (:,:,3) ) 
     118      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'      , u_ice     (:,:)   ) 
     119      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'      , v_ice     (:,:)   ) 
     120      CALL iom_rstput( iter, nitrst, numriw, 'qstoif'     , qstoif    (:,:)   ) 
     121      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'      , fsbbq     (:,:)   ) 
     122#if ! defined key_lim2_vp 
     123      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'  , stress1_i (:,:)   ) 
     124      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'  , stress2_i (:,:)   ) 
     125      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i' , stress12_i(:,:)   ) 
     126#endif 
     127      CALL iom_rstput( iter, nitrst, numriw, 'sxice' ,      sxice     (:,:)   ) 
     128      CALL iom_rstput( iter, nitrst, numriw, 'syice' ,      syice     (:,:)   ) 
     129      CALL iom_rstput( iter, nitrst, numriw, 'sxxice',      sxxice    (:,:)   ) 
     130      CALL iom_rstput( iter, nitrst, numriw, 'syyice',      syyice    (:,:)   ) 
     131      CALL iom_rstput( iter, nitrst, numriw, 'sxyice',      sxyice    (:,:)   ) 
     132      CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  ,      sxsn      (:,:)   ) 
     133      CALL iom_rstput( iter, nitrst, numriw, 'sysn'  ,      sysn      (:,:)   ) 
     134      CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' ,      sxxsn     (:,:)   ) 
     135      CALL iom_rstput( iter, nitrst, numriw, 'syysn' ,      syysn     (:,:)   ) 
     136      CALL iom_rstput( iter, nitrst, numriw, 'sxysn' ,      sxysn     (:,:)   ) 
     137      CALL iom_rstput( iter, nitrst, numriw, 'sxa'   ,      sxa       (:,:)   ) 
     138      CALL iom_rstput( iter, nitrst, numriw, 'sya'   ,      sya       (:,:)   ) 
     139      CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  ,      sxxa      (:,:)   ) 
     140      CALL iom_rstput( iter, nitrst, numriw, 'syya'  ,      syya      (:,:)   ) 
     141      CALL iom_rstput( iter, nitrst, numriw, 'sxya'  ,      sxya      (:,:)   ) 
     142      CALL iom_rstput( iter, nitrst, numriw, 'sxc0'  ,      sxc0      (:,:)   ) 
     143      CALL iom_rstput( iter, nitrst, numriw, 'syc0'  ,      syc0      (:,:)   ) 
     144      CALL iom_rstput( iter, nitrst, numriw, 'sxxc0' ,      sxxc0     (:,:)   ) 
     145      CALL iom_rstput( iter, nitrst, numriw, 'syyc0' ,      syyc0     (:,:)   ) 
     146      CALL iom_rstput( iter, nitrst, numriw, 'sxyc0' ,      sxyc0     (:,:)   ) 
     147      CALL iom_rstput( iter, nitrst, numriw, 'sxc1'  ,      sxc1      (:,:)   ) 
     148      CALL iom_rstput( iter, nitrst, numriw, 'syc1'  ,      syc1      (:,:)   ) 
     149      CALL iom_rstput( iter, nitrst, numriw, 'sxxc1' ,      sxxc1     (:,:)   ) 
     150      CALL iom_rstput( iter, nitrst, numriw, 'syyc1' ,      syyc1     (:,:)   ) 
     151      CALL iom_rstput( iter, nitrst, numriw, 'sxyc1' ,      sxyc1     (:,:)   ) 
     152      CALL iom_rstput( iter, nitrst, numriw, 'sxc2'  ,      sxc2      (:,:)   ) 
     153      CALL iom_rstput( iter, nitrst, numriw, 'syc2'  ,      syc2      (:,:)   ) 
     154      CALL iom_rstput( iter, nitrst, numriw, 'sxxc2' ,      sxxc2     (:,:)   ) 
     155      CALL iom_rstput( iter, nitrst, numriw, 'syyc2' ,      syyc2     (:,:)   ) 
     156      CALL iom_rstput( iter, nitrst, numriw, 'sxyc2' ,      sxyc2     (:,:)   ) 
     157      CALL iom_rstput( iter, nitrst, numriw, 'sxst'  ,      sxst      (:,:)   ) 
     158      CALL iom_rstput( iter, nitrst, numriw, 'syst'  ,      syst      (:,:)   ) 
     159      CALL iom_rstput( iter, nitrst, numriw, 'sxxst' ,      sxxst     (:,:)   ) 
     160      CALL iom_rstput( iter, nitrst, numriw, 'syyst' ,      syyst     (:,:)   ) 
     161      CALL iom_rstput( iter, nitrst, numriw, 'sxyst' ,      sxyst     (:,:)   ) 
    156162       
    157163      IF( iter == nitrst ) THEN 
    … …  
    218224      ENDIF 
    219225 
    220       CALL iom_get( numrir, jpdom_autoglo, 'qstoif', qstoif )     
    221       CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq  )     
    222       CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  ) 
    223       CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  ) 
    224       CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice ) 
    225       CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice ) 
    226       CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice ) 
    227       CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   ) 
    228       CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   ) 
    229       CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  ) 
    230       CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  ) 
    231       CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  ) 
    232       CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    ) 
    233       CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    ) 
    234       CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   ) 
    235       CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   ) 
    236       CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   ) 
    237       CALL iom_get( numrir, jpdom_autoglo, 'sxc0'  , sxc0   ) 
    238       CALL iom_get( numrir, jpdom_autoglo, 'syc0'  , syc0   ) 
    239       CALL iom_get( numrir, jpdom_autoglo, 'sxxc0' , sxxc0  ) 
    240       CALL iom_get( numrir, jpdom_autoglo, 'syyc0' , syyc0  ) 
    241       CALL iom_get( numrir, jpdom_autoglo, 'sxyc0' , sxyc0  ) 
    242       CALL iom_get( numrir, jpdom_autoglo, 'sxc1'  , sxc1   ) 
    243       CALL iom_get( numrir, jpdom_autoglo, 'syc1'  , syc1   ) 
    244       CALL iom_get( numrir, jpdom_autoglo, 'sxxc1' , sxxc1  ) 
    245       CALL iom_get( numrir, jpdom_autoglo, 'syyc1' , syyc1  ) 
    246       CALL iom_get( numrir, jpdom_autoglo, 'sxyc1' , sxyc1  ) 
    247       CALL iom_get( numrir, jpdom_autoglo, 'sxc2'  , sxc2   ) 
    248       CALL iom_get( numrir, jpdom_autoglo, 'syc2'  , syc2   ) 
    249       CALL iom_get( numrir, jpdom_autoglo, 'sxxc2' , sxxc2  ) 
    250       CALL iom_get( numrir, jpdom_autoglo, 'syyc2' , syyc2  ) 
    251       CALL iom_get( numrir, jpdom_autoglo, 'sxyc2' , sxyc2  ) 
    252       CALL iom_get( numrir, jpdom_autoglo, 'sxst'  , sxst   ) 
    253       CALL iom_get( numrir, jpdom_autoglo, 'syst'  , syst   ) 
    254       CALL iom_get( numrir, jpdom_autoglo, 'sxxst' , sxxst  ) 
    255       CALL iom_get( numrir, jpdom_autoglo, 'syyst' , syyst  ) 
    256       CALL iom_get( numrir, jpdom_autoglo, 'sxyst' , sxyst  ) 
     226      CALL iom_get( numrir, jpdom_autoglo, 'qstoif'     , qstoif     )     
     227      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'      , fsbbq      )     
     228#if ! defined key_lim2_vp 
     229      CALL iom_get( numrir, jpdom_autoglo, 'stress1_i'  , stress1_i  ) 
     230      CALL iom_get( numrir, jpdom_autoglo, 'stress2_i'  , stress2_i  ) 
     231      CALL iom_get( numrir, jpdom_autoglo, 'stress12_i' , stress12_i ) 
     232#endif 
     233      CALL iom_get( numrir, jpdom_autoglo, 'sxice'      , sxice      ) 
     234      CALL iom_get( numrir, jpdom_autoglo, 'syice'      , syice      ) 
     235      CALL iom_get( numrir, jpdom_autoglo, 'sxxice'     , sxxice     ) 
     236      CALL iom_get( numrir, jpdom_autoglo, 'syyice'     , syyice     ) 
     237      CALL iom_get( numrir, jpdom_autoglo, 'sxyice'     , sxyice     ) 
     238      CALL iom_get( numrir, jpdom_autoglo, 'sxsn'       , sxsn       ) 
     239      CALL iom_get( numrir, jpdom_autoglo, 'sysn'       , sysn       ) 
     240      CALL iom_get( numrir, jpdom_autoglo, 'sxxsn'      , sxxsn      ) 
     241      CALL iom_get( numrir, jpdom_autoglo, 'syysn'      , syysn      ) 
     242      CALL iom_get( numrir, jpdom_autoglo, 'sxysn'      , sxysn      ) 
     243      CALL iom_get( numrir, jpdom_autoglo, 'sxa'        , sxa        ) 
     244      CALL iom_get( numrir, jpdom_autoglo, 'sya'        , sya        ) 
     245      CALL iom_get( numrir, jpdom_autoglo, 'sxxa'       , sxxa       ) 
     246      CALL iom_get( numrir, jpdom_autoglo, 'syya'       , syya       ) 
     247      CALL iom_get( numrir, jpdom_autoglo, 'sxya'       , sxya       ) 
     248      CALL iom_get( numrir, jpdom_autoglo, 'sxc0'       , sxc0       ) 
     249      CALL iom_get( numrir, jpdom_autoglo, 'syc0'       , syc0       ) 
     250      CALL iom_get( numrir, jpdom_autoglo, 'sxxc0'      , sxxc0      ) 
     251      CALL iom_get( numrir, jpdom_autoglo, 'syyc0'      , syyc0      ) 
     252      CALL iom_get( numrir, jpdom_autoglo, 'sxyc0'      , sxyc0      ) 
     253      CALL iom_get( numrir, jpdom_autoglo, 'sxc1'       , sxc1       ) 
     254      CALL iom_get( numrir, jpdom_autoglo, 'syc1'       , syc1       ) 
     255      CALL iom_get( numrir, jpdom_autoglo, 'sxxc1'      , sxxc1      ) 
     256      CALL iom_get( numrir, jpdom_autoglo, 'syyc1'      , syyc1      ) 
     257      CALL iom_get( numrir, jpdom_autoglo, 'sxyc1'      , sxyc1      ) 
     258      CALL iom_get( numrir, jpdom_autoglo, 'sxc2'       , sxc2       ) 
     259      CALL iom_get( numrir, jpdom_autoglo, 'syc2'       , syc2       ) 
     260      CALL iom_get( numrir, jpdom_autoglo, 'sxxc2'      , sxxc2      ) 
     261      CALL iom_get( numrir, jpdom_autoglo, 'syyc2'      , syyc2      ) 
     262      CALL iom_get( numrir, jpdom_autoglo, 'sxyc2'      , sxyc2      ) 
     263      CALL iom_get( numrir, jpdom_autoglo, 'sxst'       , sxst       ) 
     264      CALL iom_get( numrir, jpdom_autoglo, 'syst'       , syst       ) 
     265      CALL iom_get( numrir, jpdom_autoglo, 'sxxst'      , sxxst      ) 
     266      CALL iom_get( numrir, jpdom_autoglo, 'syyst'      , syyst      ) 
     267      CALL iom_get( numrir, jpdom_autoglo, 'sxyst'      , sxyst      ) 
    257268       
    258269      CALL iom_close( numrir ) 

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/limsbc_2.F90

    r1756 r2135  
    44   !!           computation of the flux at the sea ice/ocean interface 
    55   !!====================================================================== 
    6    !! History : 00-01 (H. Goosse) Original code 
    7    !!           02-07 (C. Ethe, G. Madec) re-writing F90 
    8    !!           06-07 (G. Madec) surface module 
     6   !! History :  LIM  ! 2000-01 (H. Goosse) Original code 
     7   !!            1.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90 
     8   !!            3.0  ! 2006-07 (G. Madec) surface module 
     9   !!            3.3  ! 2009-05 (G.Garric, C. Bricaud) addition of the lim2_evp case 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_lim2 
    1112   !!---------------------------------------------------------------------- 
    1213   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    13    !!---------------------------------------------------------------------- 
    1414   !!---------------------------------------------------------------------- 
    1515   !!   lim_sbc_2  : flux at the ice / ocean interface 
    … …  
    1717   USE par_oce          ! ocean parameters 
    1818   USE dom_oce          ! ocean domain 
    19    USE sbc_ice          ! surface boundary condition 
    20    USE sbc_oce          ! surface boundary condition 
     19   USE sbc_ice          ! surface boundary condition: ice 
     20   USE sbc_oce          ! surface boundary condition: ocean 
    2121   USE phycst           ! physical constants 
    22    USE ice_2            ! LIM sea-ice variables 
    23  
    24    USE lbclnk           ! ocean lateral boundary condition 
     22   USE ice_2            ! LIM2: ice variables 
     23 
     24   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges 
    2525   USE in_out_manager   ! I/O manager 
    2626   USE diaar5, ONLY :   lk_diaar5 
    27    USE iom              !  
     27   USE iom              ! IOM library 
    2828   USE albedo           ! albedo parameters 
    2929   USE prtctl           ! Print control 
    … …  
    3333   PRIVATE 
    3434 
    35    PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2 
    36  
    37    REAL(wp)  ::   epsi16 = 1.e-16  ! constant values 
    38    REAL(wp)  ::   rzero  = 0.e0     
    39    REAL(wp)  ::   rone   = 1.e0 
    40    REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r 
    41    REAL(wp), DIMENSION(jpi,jpj)  ::   sice_r 
     35   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2 
     36 
     37   REAL(wp)  ::   r1_rdtice                    ! constant values 
     38   REAL(wp)  ::   epsi16 = 1.e-16              !     -      - 
     39   REAL(wp)  ::   rzero  = 0.e0                !     -      - 
     40   REAL(wp)  ::   rone   = 1.e0                !     -      - 
     41   ! 
     42   REAL(wp), DIMENSION(jpi,jpj) ::   soce_r, sice_r   ! constant SSS and ice salinity used in levitating sea-ice case 
    4243 
    4344   !! * Substitutions 
    4445#  include "vectopt_loop_substitute.h90" 
    4546   !!---------------------------------------------------------------------- 
    46    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     47   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4748   !! $Id$ 
    4849   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4950   !!---------------------------------------------------------------------- 
    50  
    5151CONTAINS 
    5252 
    … …  
    7878      !! 
    7979      INTEGER  ::   ji, jj           ! dummy loop indices 
    80       INTEGER  ::   ifvt, i1mfr, idfr               ! some switches 
    81       INTEGER  ::   iflt, ial, iadv, ifral, ifrdv 
    82       INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers 
    83       REAL(wp) ::   zrdtir           ! 1. / rdt_ice 
    84       REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux 
    85       REAL(wp) ::   zinda            ! switch for testing the values of ice concentration 
    86       REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface 
    87       REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface 
    88       REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points 
    89       REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points 
    90       REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity 
    91 ! interface 2D --> 3D 
    92       REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb     ! albedo of ice under overcast sky 
    93       REAL(wp), DIMENSION(jpi,jpj,1) ::   zalbp    ! albedo of ice under clear sky 
    94       REAL(wp) ::   zsang, zmod, zztmp, zfm 
    95       REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! component of ocean stress below sea-ice at I-point 
    96       REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi           ! module    of ocean stress below sea-ice at I-point 
    97       REAL(wp), DIMENSION(jpi,jpj) ::   zqnsoce          ! save qns before its modification by ice model 
    98  
     80      INTEGER  ::   ii0, ii1, ij0, ij1         ! local integers 
     81      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       - 
     82      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       - 
     83      REAL(wp) ::   zqsr, zqns, zsang, zmod, zfm   ! local scalars 
     84      REAL(wp) ::   zinda, zfons, zemp, zztmp      !   -      - 
     85      REAL(wp) ::   zfrldu, zutau, zu_io           !   -      - 
     86      REAL(wp) ::   zfrldv, zvtau, zv_io           !   -      - 
     87      REAL(wp), DIMENSION(jpi,jpj)   ::   ztio_u, ztio_v    ! 2D workspace 
     88      REAL(wp), DIMENSION(jpi,jpj)   ::   ztiomi, zqnsoce   !  -     - 
     89      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp   ! 2D/3D workspace 
    9990      !!--------------------------------------------------------------------- 
    10091      
    101       zrdtir = 1. / rdt_ice 
    10292       
    10393      IF( kt == nit000 ) THEN 
    … …  
    10595         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition' 
    10696         IF(lwp) WRITE(numout,*) '~~~~~~~~~   ' 
    107  
     97         ! 
     98         r1_rdtice = 1. / rdt_ice 
     99         ! 
    108100         soce_r(:,:) = soce 
    109101         sice_r(:,:) = sice 
    110102         ! 
    111          IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN 
    112             !                                        ! ======================= 
    113             !                                        !  ORCA_R2 configuration 
    114             !                                        ! ======================= 
     103         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN     !  ORCA_R2 configuration 
    115104            ii0 = 145   ;   ii1 = 180        ! Baltic Sea 
    116105            ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 
    … …  
    175164!!$ 
    176165 
    177             !   computation the solar flux at ocean surface 
     166            ! solar flux at ocean surface 
    178167#if defined key_coupled  
    179168            zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) ) 
    … …  
    181170            zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj) 
    182171#endif             
    183             !  computation the non solar heat flux at ocean surface 
     172            ! non solar heat flux at ocean surface 
    184173            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads 
    185                &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            & 
    186                &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    & 
    187                &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir 
    188  
     174               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                              & 
     175               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   & 
     176               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice 
     177            ! 
    189178            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ??? 
    190              
     179            ! 
    191180            qsr  (ji,jj) = zqsr                                          ! solar heat flux  
    192181            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux 
    … …  
    194183      END DO 
    195184 
    196       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )       
    197       CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )       
    198       CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) ) 
     185      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:)                         )       
     186      CALL iom_put( 'qns_io_cea'  , qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )       
     187      CALL iom_put( 'qsr_io_cea'  , fstric(:,:) * ( 1.e0 - pfrld(:,:) ) ) 
    199188 
    200189      !------------------------------------------! 
    201190      !      mass flux at the ocean surface      ! 
    202191      !------------------------------------------! 
    203  
    204 !!gm 
    205 !!gm CAUTION    
    206 !!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm 
    207 !!gm 
    208192      DO jj = 1, jpj 
    209193         DO ji = 1, jpi 
    210              
    211194#if defined key_coupled 
    212           zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  
    213              &   + rdmsnif(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting  
     195            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode) 
     196            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )   &   ! atmosphere-ocean freshwater flux 
     197               &                  + rdmsnif(ji,jj) * r1_rdtice                   ! freshwater flux due to snow melting  
    214198#else 
    215 !!$            !  computing freshwater exchanges at the ice/ocean interface 
    216 !!$            zpme = - evap(ji,jj)    *   frld(ji,jj)           &   !  evaporation over oceanic fraction 
    217 !!$               &   + tprecip(ji,jj)                           &   !  total precipitation 
    218 !!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice 
    219 !!$               &   - rdmsnif(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting  
    220             !  computing freshwater exchanges at the ice/ocean interface 
    221             zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction  
    222                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean 
    223                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step 
    224                &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting  
    225                !                                                   !  ice-covered fraction: 
     199            ! freshwater exchanges at the ice-atmosphere / ocean interface (forced mode) 
     200            zemp = + emp(ji,jj)     *         frld(ji,jj)     &   ! e-p budget over open ocean fraction  
     201               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )   &   ! liquid precipitation reaches directly the ocean 
     202               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   ! (account for change in ice cover within the timestep 
     203               &   + rdmsnif(ji,jj) * r1_rdtice                   ! freshwaterflux due to snow melting  
    226204#endif             
    227  
    228             !  computing salt exchanges at the ice/ocean interface 
    229             zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir )  
    230              
    231             !  converting the salt flux from ice to a freshwater flux from ocean 
     205            ! salt exchanges at the ice/ocean interface 
     206            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice )  
     207            ! 
     208            ! convert the salt flux from ice into a freshwater flux from ocean 
    232209            zfm  = zfons / ( sss_m(ji,jj) + epsi16 ) 
    233              
     210            ! 
    234211            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution) 
    235212            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution) 
    236  
    237213         END DO 
    238214      END DO 
    239  
     215      ! 
    240216      IF( lk_diaar5 ) THEN 
    241          CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * zrdtir ) 
    242          CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * zrdtir ) 
    243          CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * zrdtir ) 
     217         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice ) 
     218         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * r1_rdtice ) 
     219         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * r1_rdtice ) 
    244220      ENDIF 
    245221 
    … …  
    275251            DO ji = 2, jpim1   ! NO vector opt. 
    276252               ! ... components of ice-ocean stress at U and V-points  (from I-point values) 
    277                zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) ) 
     253#if defined key_lim2_vp 
     254               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )      ! VP rheology 
    278255               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) ) 
     256#else 
     257               zutau  = ztio_u(ji,jj)                                      ! EVP rheology 
     258               zvtau  = ztio_v(ji,jj) 
     259#endif 
    279260               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) 
    280261               zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj  ) ) 
    … …  
    290271            END DO 
    291272         END DO 
    292  
    293          ! boundary condition on the stress (utau,vtau,taum) 
    294          CALL lbc_lnk( utau, 'U', -1. ) 
    295          CALL lbc_lnk( vtau, 'V', -1. ) 
     273         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition  
    296274         CALL lbc_lnk( taum, 'T',  1. ) 
    297275 
    298276      ENDIF 
    299277 
     278      IF( lk_cpl ) THEN            
    300279      !-----------------------------------------------! 
    301280      !   Coupling variables                          ! 
    302281      !-----------------------------------------------! 
    303  
    304       IF ( lk_cpl ) THEN            
    305          ! Ice surface temperature  
    306          tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature        
    307          ! Computation of snow/ice and ocean albedo 
     282         tn_ice(:,:,1) = sist(:,:)           ! sea-ice surface temperature        
     283         !                                   ! snow/ice and ocean albedo 
    308284         CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb ) 
    309285         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys) 
    … …  
    318294         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ') 
    319295      ENDIF  
    320     
    321     END SUBROUTINE lim_sbc_2 
     296      ! 
     297   END SUBROUTINE lim_sbc_2 
    322298 
    323299#else 

  • branches/devmercator2010_1/NEMO/LIM_SRC_2/limtrp_2.F90

    r1715 r2135  
    44   !! LIM 2.0 transport ice model : sea-ice advection/diffusion 
    55   !!====================================================================== 
     6   !! History :  LIM  !  2000-01  (LIM)  Original code 
     7   !!            2.0  !  2001-05  (G. Madec, R. Hordoir) opa norm 
     8   !!             -   !  2004-01  (G. Madec, C. Ethe)  F90, mpp 
     9   !!            3.3  !  2009-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
     10   !!---------------------------------------------------------------------- 
    611#if defined key_lim2 
    712   !!---------------------------------------------------------------------- 
    … …  
    1116   !!   lim_trp_init_2 : initialization and namelist read 
    1217   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    14    USE phycst 
    15    USE dom_oce 
     18   USE phycst          ! physical constants 
     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       ! LIM2 sea-ice domain 
     22   USE ice_2           ! LIM2 variables 
     23   USE limistate_2     ! LIM2 initial state 
     24   USE limadv_2        ! LIM2 advection 
     25   USE limhdf_2        ! LIM2 horizontal diffusion 
     26   USE lbclnk          ! Lateral Boundary condition - 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 
    36    REAL(wp)  ::           &  ! constant values 
    37       epsi06 = 1.e-06  ,  & 
    38       epsi03 = 1.e-03  ,  & 
    39       epsi16 = 1.e-16  ,  & 
    40       rzero  = 0.e0    ,  & 
    41       rone   = 1.e0 
     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 
     36   REAL(wp)  ::   epsi06 = 1.e-06   ! constant values 
     37   REAL(wp)  ::   epsi03 = 1.e-03   
     38   REAL(wp)  ::   epsi16 = 1.e-16   
     39   REAL(wp)  ::   rzero  = 0.e0    
     40   REAL(wp)  ::   rone   = 1.e0 
    4241 
    4342   !! * Substitution 
    4443#  include "vectopt_loop_substitute.h90" 
    4544   !!---------------------------------------------------------------------- 
    46    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
     45   !! NEMO/LIM2 3.3,  UCL-LOCEAN-IPSL (2010)  
    4746   !! $Id$ 
    48    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    49    !!---------------------------------------------------------------------- 
    50  
     47   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     48   !!---------------------------------------------------------------------- 
    5149CONTAINS 
    5250 
    … …  
    6260      !! 
    6361      !! ** 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 
    6962      !!--------------------------------------------------------------------- 
    7063      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 
     64      !! 
     65      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     66      INTEGER  ::   initad       ! number of sub-timestep for the advection 
     67 
     68      REAL(wp) ::   zindb  , zacrith, zindsn , zindic , zusvosn   ! local scalars 
     69      REAL(wp) ::   zusvoic, zignm  , zindhe , zvbord , zcfl      !   -      - 
     70      REAL(wp) ::   zusnit , zrtt   , ztsn   , ztic1  , ztic2     !   -      - 
     71 
     72      REAL(wp), DIMENSION(jpi,jpj) ::   zui_u , zvi_v , zsm             ! 2D workspace 
     73      REAL(wp), DIMENSION(jpi,jpj) ::   zs0ice, zs0sn , zs0a            !  -      - 
     74      REAL(wp), DIMENSION(jpi,jpj) ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      - 
    9475      !--------------------------------------------------------------------- 
    9576 
    … …  
    10788         ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.         
    10889         zvbord = 1.0 + ( 1.0 - bound ) 
    109          DO jj = 1, jpjm1 
     90#if defined key_lim2_vp 
     91         DO jj = 1, jpjm1     ! VP rheology: ice (u,v) at I-point 
    11092            DO ji = 1, jpim1   ! NO vector opt. 
    111                zui_u(ji,jj) = ( u_ice(ji+1,jj  ) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) ) 
    112                zvi_v(ji,jj) = ( v_ice(ji  ,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) ) 
     93               zui_u(ji,jj) = ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) ) 
     94               zvi_v(ji,jj) = ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) ) 
    11395            END DO 
    11496         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. ) 
     97         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )      ! Lateral boundary conditions  
     98#else 
     99        zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points 
     100        zvi_v(:,:) = v_ice(:,:) 
     101#endif 
    118102 
    119103         ! CFL test for stability 
    … …  
    123107         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) ) 
    124108 
    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 
     109         IF(lk_mpp ) CALL mpp_max(zcfl) 
     110 
     111         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : cfl criterion violation. day ',nday,' cfl = ',zcfl 
    128112 
    129113         ! content of properties 
    130114         ! --------------------- 
    131          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. 
    135          zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer. 
    136          zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer. 
    137          zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a(:,:)    ! Heat reservoir for brine pockets. 
    138           
     115         zs0sn (:,:) =  hsnm(:,:)              * area(:,:)     ! Snow volume. 
     116         zs0ice(:,:) =  hicm (:,:)             * area(:,:)     ! Ice volume. 
     117         zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area(:,:)     ! Surface covered by ice. 
     118         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn(:,:)    ! Heat content of the snow layer. 
     119         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)   ! Heat content of the first ice layer. 
     120         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)   ! Heat content of the second ice layer. 
     121         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)   ! Heat reservoir for brine pockets. 
    139122  
    140123         ! Advection  
    … …  
    144127         zusnit = 1.0 / REAL( initad )  
    145128          
    146          IF ( MOD( nday , 2 ) == 0) THEN 
     129!!gm this has been changed in the reference to become odd and even ice time step 
     130 
     131         IF( MOD( nday , 2 ) == 0) THEN 
    147132            DO jk = 1,initad 
    148133               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) 
    … …  
    231216         ! Up-dating and limitation of sea ice properties after transport. 
    232217         DO jj = 1, jpj 
    233 !!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq! 
    234218            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH 
    235219            DO ji = 1, jpi 
    … …  
    272256            END DO 
    273257         END DO 
    274           
     258         ! 
    275259      ENDIF 
    276        
     260      ! 
    277261   END SUBROUTINE lim_trp_2 
    278262 
    … …  
    288272      !! 
    289273      !! ** input   :   Namelist namicetrp 
    290       !! 
    291       !! history : 
    292       !!   2.0  !  03-08 (C. Ethe)  Original code 
    293274      !!------------------------------------------------------------------- 
    294275      NAMELIST/namicetrp/ bound 
    295276      !!------------------------------------------------------------------- 
    296  
     277      ! 
    297278      ! Read Namelist namicetrp 
    298279      REWIND ( numnam_ice ) 
    … …  
    304285         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound 
    305286      ENDIF 
    306              
     287      !      
    307288   END SUBROUTINE lim_trp_init_2 
    308289 
Note: See TracChangeset for help on using the changeset viewer.