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 2095 – NEMO

Changeset 2095


Ignore:
Timestamp:
2010-09-15T14:10:33+02:00 (14 years ago)
Author:
cbricaud
Message:

add correction from reviewer

Location:
branches/dev_1784_EVP/NEMO
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • branches/dev_1784_EVP/NEMO/LIM_SRC_2/dom_ice_2.F90

    r2046 r2095  
    44   !! LIM 2.0 Sea Ice :   Domain  variables 
    55   !!====================================================================== 
    6    !! History :   2.0  !  03-08  (C. Ethe)  Free form and module 
    7    !!             3.3  !  09-05 (G.Garric) addition of the lim2_evp cas 
     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 
    88   !!---------------------------------------------------------------------- 
    9 #if defined key_lim2 
     9#if defined   key_lim2 
    1010   !!---------------------------------------------------------------------- 
    11    !!   LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    12    !! $Id$ 
    13    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     11   !!   'key_lim2'                                       LIM2 sea-ice model 
    1412   !!---------------------------------------------------------------------- 
    1513   USE par_ice_2 
     
    2119 
    2220   INTEGER, PUBLIC ::   njeq , njeqm1          !: j-index of the equator if it is inside the domain 
    23       !                                        !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
     21   !                                           !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
    2422 
    25    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor,   &  !: coriolis factor and coeficient 
    26       &                                              covrai ,         &  !: sine of geographic latitude 
    27       &                                              area   ,         &  !: surface of grid cell  
    28       &                                              tms    , tmu        !: temperature and velocity points masks 
    29    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght                !: weight of the 4 neighbours to compute averages 
     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 
    3028 
    31 #if defined key_lim2_vp 
    32    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::                     & 
    33       &                                              akappa , bkappa     !: first and third group of metric coefficients 
    34    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd              !: second group of metric coefficients 
     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 
    3537#else 
    36    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmv    , tmf        !: y-velocity and F-points masks 
    37    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmi                 !: ice mask: =1 if ice thick > 0 
     38   !!---------------------------------------------------------------------- 
     39   !!   Default option          Empty module         NO LIM-2 sea-ice model 
     40   !!---------------------------------------------------------------------- 
    3841#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) 
    3947   !!====================================================================== 
    40 #endif 
    4148END MODULE dom_ice_2 
  • branches/dev_1784_EVP/NEMO/LIM_SRC_2/ice_2.F90

    r2046 r2095  
    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 
    7    !!            3.3  !  09-05 (G.Garric) addition of the lim2_evp cas 
     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 
    88   !!---------------------------------------------------------------------- 
    99#if defined key_lim2 
     
    1111   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
    13    !!  LIM 2.0, UCL-LOCEAN-IPSL (2006) 
    14    !! $Id$ 
    15    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    16    !!---------------------------------------------------------------------- 
    17    !! * Modules used 
    1813   USE par_ice_2          ! LIM sea-ice parameters 
    1914 
     
    2116   PRIVATE 
    2217 
    23    !!* Share parameters namelist (namicerun read in iceini) * 
     18   !                                                                     !!* namicerun namelist * 
    2419   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input) 
    2520   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output) 
    2621   LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F) 
    2722   LOGICAL               , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping 
    28    LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag for sea-ice points output (T) or not (F) 
    29    REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of temp. in snow (0) or not (9999) 
    30    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) 
    3126   REAL(wp), DIMENSION(2), PUBLIC ::   acrit = (/ 1.e-06 , 1.e-06 /)      !: minimum fraction for leads in  
    3227   !                                                                      !: north and south hemisphere 
    33    !!* ice-dynamic namelist (namicedyn) * 
     28   !                                       !!* ice-dynamic namelist (namicedyn) * 
    3429   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation 
    3530   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation 
     
    5651   REAL(wp), PUBLIC ::   pstarh             !: pstar / 2.0 
    5752 
    58    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv         !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
    59    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv         !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
    60    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s               !: friction velocity 
     53   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
     54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
     55   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s         !: friction velocity 
    6156 
    62 #if  defined key_lim2_vp 
    63    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm         !: mean snow and ice thicknesses 
     57# if defined key_lim2_vp 
     58   !                                                        !!* VP rheology * 
     59   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm     !: mean snow and ice thicknesses 
    6460   CHARACTER(len=1), PUBLIC :: cl_grid =   'B'                   !: type of grid used in ice dynamics, 'C' or 'B' 
    65 #else 
    66    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::                      & 
    67                                                    stress1_i ,  &!: first stress tensor element        
    68                                                    stress2_i ,  &!: second stress tensor element 
    69                                                    stress12_i,  &!: diagonal stress tensor element 
    70                                                    delta_i   ,  &!: Delta factor for the ice rheology (see Flato and Hibler 95) [s-1] -> limrhg.F90 
    71                                                    divu_i    ,  &!: Divergence of the velocity field [s-1] -> limrhg.F90 
    72                                                    shear_i       !: Shear of the velocity field [s-1] -> limrhg.F90 
    73  
    74    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         :: at_i          !: 
    75    REAL(wp), PUBLIC, DIMENSION(:,:)    ,POINTER :: vt_s ,vt_i    !: mean snow and ice thicknesses 
    76    REAL(wp), PUBLIC, DIMENSION(jpi,jpj),TARGET  :: hsnm , hicm   !: mean snow and ice thicknesses, target for pointers vt_s and vt_i  
    77    CHARACTER(len=1), PUBLIC :: cl_grid =   'C'                   !: 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  
    7874#endif 
    7975 
    80    !!* diagnostic quantities 
    81    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvosif       !: Variation of volume at surface (only used for outputs) 
    82    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvobif       !: Variation of ice volume at the bottom ice (only used for outputs) 
    83    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdvolif       !: Total variation of ice volume (only used for outputs) 
    84    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvonif       !: Lateral Variation of ice volume (only used for outputs) 
    85    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin) 
     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] 
    8682   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS 
    8783   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness 
     
    9692   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass 
    9793   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean) 
    98    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  
    9995   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean 
    10096   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2) 
    10197   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget 
    10298   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice 
    103    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 
    104100   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?) 
    105101   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?) 
    106    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) 
    107103   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice 
    108104 
    109    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice   !: two components of the ice   velocity at I-point (m/s) 
    110    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] 
    111107 
    112108   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer 
    113109 
    114    !!* moment used in the advection scheme 
     110   !                                                                               !!* moment used in the advection 
    115111   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume 
    116112   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume 
     
    123119#else 
    124120   !!---------------------------------------------------------------------- 
    125    !!   Default option         Empty module        NO LIM 2.0 sea-ice model 
     121   !!   Default option          Empty module         NO LIM-2 sea-ice model 
    126122   !!---------------------------------------------------------------------- 
    127123#endif 
    128124 
     125   !!---------------------------------------------------------------------- 
     126   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
     127   !! $Id$ 
     128   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    129129   !!====================================================================== 
    130130END MODULE ice_2 
  • branches/dev_1784_EVP/NEMO/LIM_SRC_2/iceini_2.F90

    r2046 r2095  
    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) addition of the lim2_evp case 
     8   !!             3.3  !  09-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
    99   !!---------------------------------------------------------------------- 
    1010#if defined key_lim2 
     
    1212   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!---------------------------------------------------------------------- 
    1514   !!   ice_init_2       : sea-ice model initialization 
    1615   !!   ice_run_2        : Definition some run parameter for ice model 
    1716   !!---------------------------------------------------------------------- 
    18    USE dom_oce 
    19    USE dom_ice_2 
    20    USE sbc_oce         ! surface boundary condition: ocean 
    21    USE sbc_ice         ! surface boundary condition: ice 
    22    USE phycst          ! Define parameters for the routines 
    23    USE ice_2 
    24    USE limmsh_2 
    25    USE limistate_2 
    26    USE limrst_2    
    27    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 
    2827       
    2928   IMPLICIT NONE 
    3029   PRIVATE 
    3130 
    32    PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90 
     31   PUBLIC   ice_init_2   ! called by sbcice_lim_2.F90 
    3332 
    34    INTEGER , PUBLIC  ::     numit                  !: iteration number 
    35  
     33   INTEGER, PUBLIC ::   numit   !: iteration number 
    3634 
    3735   !!---------------------------------------------------------------------- 
    38    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    39    !! $Id$  
    40    !! 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) 
    4139   !!---------------------------------------------------------------------- 
    42  
    4340CONTAINS 
    4441 
     
    4744      !!                  ***  ROUTINE ice_init_2  *** 
    4845      !! 
    49       !! ** purpose :    
     46      !! ** purpose :   initialisation of LIM-2 domain and variables  
    5047      !!---------------------------------------------------------------------- 
    5148      ! 
    52       ! Open the namelist file  
     49      !                               ! Open the namelist file  
    5350      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )       
    54       CALL ice_run_2                    !  read in namelist some run parameters 
     51      CALL ice_run_2                  !  read in namelist some run parameters 
    5552                  
    56       ! Louvain la Neuve Ice model 
    57       rdt_ice = nn_fsbc * rdttra(1) 
     53      rdt_ice = nn_fsbc * rdttra(1)   ! sea-ice time step 
     54      numit   = nit000 - 1               
    5855 
    5956      CALL lim_msh_2                  ! ice mesh initialization 
    6057      
    61       ! Initial sea-ice state 
    62       IF( .NOT.ln_rstart ) THEN 
    63          CALL lim_istate_2            ! start from rest: sea-ice deduced from sst 
    64       ELSE 
    65          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 
    6661      ENDIF 
    6762       
    6863      tn_ice(:,:,1) = sist(:,:)       ! initialisation of ice temperature    
    6964      fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction     
    70       ! 
    71       numit = nit000 - 1              !initialisation ice time-step 
    72  
    7365      ! 
    7466   END SUBROUTINE ice_init_2 
     
    8274      !! 
    8375      !! ** Method  :   Read the namicerun namelist and check the parameter  
    84       !!       values called at the first timestep (nit000) 
     76      !!              values called at the first timestep (nit000) 
    8577      !! 
    8678      !! ** input   :   Namelist namicerun 
     
    8981      !!------------------------------------------------------------------- 
    9082      !                     
    91       REWIND ( numnam_ice )                       ! Read Namelist namicerun  
    92       READ   ( numnam_ice , namicerun ) 
    93  
    94       IF(lwp) THEN 
     83      REWIND( numnam_ice )                   ! Read namicerun namelist 
     84      READ  ( numnam_ice , namicerun ) 
     85      ! 
     86      IF(lwp) THEN                           ! control print 
    9587         WRITE(numout,*) 
    9688         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
  • branches/dev_1784_EVP/NEMO/LIM_SRC_2/limdyn_2.F90

    r2046 r2095  
    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) addition of the lim2_evp cas 
     10   !!             3.3  !  09-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
    1111   !!--------------------------------------------------------------------- 
    1212#if defined key_lim2 
     
    1717   !!    lim_dyn_init_2 : initialization and namelist read 
    1818   !!---------------------------------------------------------------------- 
    19    USE dom_oce        ! ocean space and time domain 
    20    USE sbc_oce        ! 
    21    USE phycst         ! 
    22    USE ice_2          ! 
    23    USE dom_ice_2      ! 
    24    USE limistate_2    ! 
    25 #if defined key_lim2_vp 
    26    USE limrhg_2       ! ice rheology 
    27 #else 
    28    USE limrhg         ! ice rheology 
    29 #endif 
    30    USE lbclnk         ! 
    31    USE lib_mpp        ! 
    32    USE in_out_manager ! I/O manager 
    33    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 
    3434 
    3535   IMPLICIT NONE 
    3636   PRIVATE 
    3737 
    38    PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim 
    39  
    40    !! * Module variables 
    41    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 
    4241 
    4342#  include "vectopt_loop_substitute.h90" 
    4443   !!---------------------------------------------------------------------- 
    45    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     44   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4645   !! $Id$ 
    4746   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4847   !!---------------------------------------------------------------------- 
    49  
    5048CONTAINS 
    5149 
     
    7573      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only) 
    7674       
    77       IF( ln_limdyn ) THEN 
    78          ! 
    79          ! Mean ice and snow thicknesses.           
    80          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 
    8179         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:) 
    8280         ! 
    83          !                                     ! Rheology (ice dynamics) 
    84          !                                     ! ======== 
     81         !                                     ! Define the j-limits where ice rheology is computed 
    8582          
    86          !  Define the j-limits where ice rheology is computed 
    87          ! --------------------------------------------------- 
    88           
    89          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  ==! 
    9084            i_j1 = 1    
    9185            i_jpj = jpj 
    9286            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    9387#if defined key_lim2_vp 
    94             CALL lim_rhg_2( i_j1, i_jpj ) 
    95 #else 
    96             CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt ) 
    97 #endif 
     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 
    9897            ! 
    99          ELSE                                 ! optimization of the computational area 
    100             ! 
    101             DO jj = 1, jpj 
    102                zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line 
    103                zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line 
    104             END DO 
    105             ! 
    106             IF( l_jeq ) THEN                     ! local domain include both hemisphere 
    107                !                                 ! Rheology is computed in each hemisphere 
    108                !                                 ! only over the ice cover latitude strip 
    109                ! Northern hemisphere 
    110                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 
    111101               i_jpj = jpj 
    112102               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
     
    123113               CALL lim_rhg( i_j1, i_jpj ) 
    124114#endif 
    125                !  
    126                ! Southern hemisphere 
    127                i_j1  =  1  
     115               i_j1  =  1                    ! Southern hemisphere 
    128116               i_jpj = njeq 
    129117               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
     
    141129#endif 
    142130               !  
    143             ELSE                                 ! local domain extends over one hemisphere only 
    144                !                                 ! Rheology is computed only over the ice cover 
    145                !                                 ! latitude strip 
     131            ELSE                                      ! local domain extends over one hemisphere only: rheology is 
     132               !                                      ! computed only over the ice cover latitude strip 
    146133               i_j1  = 1 
    147134               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
     
    149136               END DO 
    150137               i_j1 = MAX( 1, i_j1-1 ) 
    151      
    152138               i_jpj  = jpj 
    153139               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    154140                  i_jpj = i_jpj - 1 
    155141               END DO 
     142               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    156143#if defined key_lim2_vp 
    157144               i_jpj = MIN( jpj, i_jpj+2) 
    158  
    159                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    160                !  
    161                CALL lim_rhg_2( i_j1, i_jpj ) 
     145               CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology 
    162146#else 
    163147               i_j1 = MAX( 1, i_j1-2 ) 
    164148               i_jpj = MIN( jpj, i_jpj+1) 
    165                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    166                CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt ) 
     149               CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology 
    167150#endif 
    168151               ! 
     
    170153            ! 
    171154         ENDIF 
    172  
     155         ! 
    173156         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 
     157         !  
    174158          
    175          ! computation of friction velocity 
    176          ! -------------------------------- 
    177  
    178      SELECT CASE( cl_grid ) 
    179  
    180       CASE( 'C' )                          ! C-grid ice dynamics 
    181          !????????????????????????????????? 
    182          ! ice-ocean velocity at U & V-points (u_ice vi_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points) 
    183          zu_io(:,:) = u_ice(:,:) - ssu_m(:,:) 
    184          zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 
    185  
    186  
    187       CASE( 'B' )                          ! B-grid ice dynamics 
    188          ! ice-ocean velocity at U & V-points (u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points) 
    189           
    190          DO jj = 1, jpjm1 
    191             DO ji = 1, jpim1   ! NO vector opt. 
    192                zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
    193                zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
    194             END DO 
    195          END DO 
    196  
    197       END SELECT 
    198  
    199          ! frictional velocity at T-point 
    200          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 
    201176            DO ji = 2, jpim1   ! NO vector opt. because of zu_io 
    202177               ust2s(ji,jj) = 0.5 * cw                                                          & 
     
    206181         END DO 
    207182         ! 
    208       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
    209          ! 
     183      ELSE                                     ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
     184         !                                     ! =============== 
    210185         zcoef = SQRT( 0.5 ) / rau0 
    211186         DO jj = 2, jpjm1 
     
    218193      ENDIF 
    219194      ! 
    220       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point 
     195      CALL lbc_lnk( ust2s, 'T',  1. )   ! lateral boundary condition 
    221196      ! 
    222197      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    223  
     198      ! 
    224199   END SUBROUTINE lim_dyn_2 
    225200 
     
    229204      !!                  ***  ROUTINE lim_dyn_init_2  *** 
    230205      !! 
    231       !! ** Purpose :   Physical constants and parameters linked to the ice 
    232       !!              dynamics 
    233       !! 
    234       !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic 
    235       !!              parameter values 
     206      !! ** Purpose :   initialisation of  the ice dynamics variables 
     207      !!               
     208      !! ** Method  :   Read the namicedyn namelist and check their values 
    236209      !! 
    237210      !! ** input   :   Namelist namicedyn 
    238211      !!------------------------------------------------------------------- 
    239       NAMELIST/namicedyn/ epsd, alpha,     & 
    240          &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
    241          &                c_rhg, etamn, creepl, ecc, ahi0, & 
    242          &                nevp, telast, alphaevp 
    243       !!------------------------------------------------------------------- 
    244  
     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      ! 
    245216      REWIND ( numnam_ice )                       ! Read Namelist namicedyn 
    246217      READ   ( numnam_ice  , namicedyn ) 
    247  
     218      ! 
    248219      IF(lwp) THEN                                ! Control print 
    249220         WRITE(numout,*) 
     
    269240         WRITE(numout,*) '       coefficient for the solution of int. stresses  alphaevp = ', alphaevp 
    270241      ENDIF 
    271  
     242      ! 
    272243      !  Initialization 
    273244      usecc2 = 1.0 / ( ecc * ecc ) 
     
    285256#else 
    286257   !!---------------------------------------------------------------------- 
    287    !!   Default option          Empty module       NO LIM 2.0 sea-ice model 
     258   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    288259   !!---------------------------------------------------------------------- 
    289260CONTAINS 
    290    SUBROUTINE lim_dyn_2         ! Empty routine 
     261   SUBROUTINE lim_dyn_2         ! Dummy routine 
    291262   END SUBROUTINE lim_dyn_2 
    292263#endif  
  • branches/dev_1784_EVP/NEMO/LIM_SRC_2/limmsh_2.F90

    r2046 r2095  
    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) 
    49       !!         additions   : 2009-05 (addition of the lim2_evp case, G. Garric) 
     45      !! References : Deleersnijder et al. Ocean Modelling 100, 7-10  
    5046      !!---------------------------------------------------------------------  
    51       !! * Local variables 
    5247      INTEGER :: ji, jj      ! dummy loop indices 
    53  
    54       REAL(wp) ::         & 
    55          zusden              ! temporary scalars 
    56 #if defined key_lim2_vp 
    57       REAL(wp), DIMENSION(jpi,jpj) ::  & 
    58          zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction 
    59          !                   ! (resp. y direction) (defined at the center) 
    60       REAL(wp) ::         & 
    61          zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid 
    62          zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid 
    63          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 
    6454#endif 
    6555      !!--------------------------------------------------------------------- 
    66  
     56      ! 
    6757      IF(lwp) THEN 
    6858         WRITE(numout,*) 
     
    7868      njeqm1 = njeq - 1  
    7969 
    80       fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor 
     70      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor at T-point 
    8171  
    82 !i    DO jj = 1, jpj 
    83 !i       zmsk(jj) = SUM( tmask(:,jj,:) )   ! = 0          if land  everywhere on a j-line 
    84 !!ii     write(numout,*) jj, zind(jj) 
    85 !i    END DO 
    86  
    8772      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere 
    8873         l_jeq = .TRUE. 
     
    129114       
    130115#if defined key_lim2_vp 
    131       ! metric coefficients for sea ice dynamic 
     116      ! metric coefficients for sea ice dynamic (VP rheology) 
    132117      !---------------------------------------- 
    133118      !                                                       ! akappa 
     
    135120         zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1) 
    136121      END DO 
    137       CALL lbc_lnk( zd1d2, 'T', -1. ) 
    138  
    139122      DO ji = 2, jpi 
    140123         zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:) 
    141124      END DO 
    142       CALL lbc_lnk( zd2d1, 'T', -1. ) 
    143  
     125      CALL lbc_lnk( zd1d2, 'T', -1. )   ;   CALL lbc_lnk( zd2d1, 'T', -1. )      ! lateral boundary condition 
     126      ! 
    144127      akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) ) 
    145128      akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) ) 
     
    147130      akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) ) 
    148131       
    149       !                                                      ! weights (wght) 
    150       DO jj = 2, jpj 
     132      !              
     133      DO jj = 2, jpj                                         ! weights (wght) at I-points 
    151134         DO ji = 2, jpi 
    152135            zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   & 
     
    163146      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. ) 
    164147#else 
    165       !                                                      ! weights (wght) 
    166       DO jj = 2, jpj-1 
    167          DO ji = 2, jpi-1 
    168             zusden = 1. / (  ( e1t(ji+1,jj) + e1t(ji,jj  ) )   & 
    169                &           * ( e2t(ji,jj+1) + e2t(ji  ,jj) ) )  
     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) ) )  
    170154            wght(ji,jj,1,1) = zusden * e1t(ji+1,jj) * e2t(ji,jj+1) 
    171155            wght(ji,jj,1,2) = zusden * e1t(ji+1,jj) * e2t(ji,jj  ) 
     
    174158         END DO 
    175159      END DO 
    176  
    177       !With EVP, the weights are calculated on 'F' points 
    178       CALL lbc_lnk( wght(:,:,1,1), 'F', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V-point 
    179       CALL lbc_lnk( wght(:,:,1,2), 'F', 1. )      ! the value of wght at jpj is wrong 
    180       CALL lbc_lnk( wght(:,:,2,1), 'F', 1. )      ! but it is never used 
    181       CALL lbc_lnk( wght(:,:,2,2), 'F', 1. ) 
    182  
     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. ) 
    183162#endif 
    184163     
    185164      ! Coefficients for divergence of the stress tensor 
    186165      !------------------------------------------------- 
    187  
    188 #if defined key_lim2_vp 
    189       DO jj = 2, jpj 
     166#if defined key_lim2_vp 
     167      DO jj = 2, jpj                                     ! VP rheology 
    190168         DO ji = 2, jpi   ! NO vector opt. 
    191             zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2)   & 
    192                &   + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    193                &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1)   & 
    194                &   + e1t(ji-1,jj-1) * wght(ji,jj,1,1) 
    195  
    196             zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2)   & 
    197                &   + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    198                &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1)   & 
    199                &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1) 
    200  
    201 ! better written but change the last digit and thus solver in less than 100 timestep 
    202 !           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   & 
    203 !              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)  
    204  
    205 !           zh2p  = e2t(ji-1,jj  ) * wght(ji,jj,1,2) + e2t(ji,jj  ) * wght(ji,jj,2,2)   & 
    206 !              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1) 
    207  
    208 !!ibug =0   zusden = 1.0 / ( zh1p * zh2p * 4.e0 ) 
    209             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 ) 
    210175            zusden2 = zusden * 2.0  
    211  
    212             zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   ) 
    213             zd2d1p = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   ) 
    214  
     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            ! 
    215179            alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1) 
    216180            alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  ) 
    217181            alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1) 
    218182            alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  ) 
    219  
     183            ! 
    220184            alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1) 
    221185            alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  ) 
    222186            alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1) 
    223187            alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  ) 
    224  
     188            ! 
    225189            alambd(ji,jj,1,2,2,1) = zd1d2p 
    226190            alambd(ji,jj,1,2,2,2) = zd1d2p 
    227191            alambd(ji,jj,1,2,1,1) = zd1d2p 
    228192            alambd(ji,jj,1,2,1,2) = zd1d2p 
    229  
     193            ! 
    230194            alambd(ji,jj,2,1,2,1) = zd2d1p 
    231195            alambd(ji,jj,2,1,2,2) = zd2d1p 
     
    234198         END DO 
    235199      END DO 
    236  
    237       CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V point 
    238       CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )      ! the value of wght at jpj is wrong 
    239       CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )      ! but it is never used 
    240       CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )      !  
    241  
    242       CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )      ! CAUTION: idem 
    243       CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )      !  
    244       CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )      ! 
    245       CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )      ! 
    246  
    247       CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )      ! CAUTION: idem 
    248       CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )      ! 
    249       CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )      ! 
    250       CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )      ! 
    251  
    252       CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )      ! CAUTION: idem 
    253       CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )      ! 
    254       CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      ! 
    255       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. ) 
    256214#endif 
    257215             
     
    259217      ! Initialization of ice masks 
    260218      !---------------------------- 
    261        
     219      !       
    262220      tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask 
    263  
    264 !i here we can use umask with a i and j shift of -1,-1 
    265       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              !  
    266225      tmu(1,:) = 0.e0 
    267  
    268 #if defined key_lim2_vp 
    269       DO jj = 2, jpj               ! ice U.V-point: computed from ice T-point mask 
     226      DO jj = 2, jpj 
    270227         DO ji = 2, jpim1   ! NO vector opt. 
    271228            tmu(ji,jj) =  tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)             
    272229         END DO 
    273230      END DO 
    274        
    275       !--lateral boundary conditions     
    276       CALL lbc_lnk( tmu(:,:), 'I', 1. ) 
    277 #else 
    278       tmv(:,1) = 0.e0 !SB 
    279       tmv(1,:) = 0.e0 !SB 
    280       tmf(1,:) = 0.e0 
    281       tmf(:,1) = 0.e0 
    282       DO jj = 1, jpj - 1 
    283          DO ji = 1 , jpi - 1 
    284             tmu(ji,jj) =  tms(ji,jj) * tms(ji+1,jj) 
    285             tmv(ji,jj) =  tms(ji,jj) * tms(ji,jj+1) 
    286             tmf(ji,jj) =  tms(ji,jj) * tms(ji+1,jj) * tms(ji,jj+1) * & 
    287                tms(ji+1,jj+1) 
    288          END DO 
    289       END DO 
    290  
    291       !--lateral boundary conditions 
    292       CALL lbc_lnk( tmu(:,:), 'U', 1. ) 
    293       CALL lbc_lnk( tmv(:,:), 'V', 1. ) 
    294       CALL lbc_lnk( tmf(:,:), 'F', 1. ) 
    295 #endif 
    296        
    297       ! unmasked and masked area of T-grid cell 
    298       area(:,:) = e1t(:,:) * e2t(:,:) 
    299        
     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      ! 
    300242   END SUBROUTINE lim_msh_2 
    301243 
  • branches/dev_1784_EVP/NEMO/LIM_SRC_2/limrhg_2.F90

    r2046 r2095  
    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    !!            3.3  !  09-05  (G.Garric) addition of the lim2_evp case 
    12    !!---------------------------------------------------------------------- 
    13 #if defined key_lim2 &&  defined key_lim2_vp 
    14    !!---------------------------------------------------------------------- 
    15    !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    16    !!---------------------------------------------------------------------- 
     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 
    1717   !!---------------------------------------------------------------------- 
    1818   !!   lim_rhg_2   : computes ice velocities 
     
    2222   USE sbc_oce        ! surface boundary condition: ocean variables 
    2323   USE sbc_ice        ! surface boundary condition: ice variables 
    24    USE dom_ice_2      ! domaine: ice variables 
    25    USE phycst         ! physical constant 
    26    USE ice_2          ! ice variables 
    27    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 
    2828   USE lib_mpp        ! MPP library 
    2929   USE in_out_manager ! I/O manager 
     
    3333   PRIVATE 
    3434 
    35    PUBLIC   lim_rhg_2 ! routine called by lim_dyn 
     35   PUBLIC   lim_rhg_2   ! routine called by lim_dyn 
    3636 
    3737   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero 
     
    4141#  include "vectopt_loop_substitute.h90" 
    4242   !!---------------------------------------------------------------------- 
    43    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     43   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4444   !! $Id$ 
    4545   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    9090      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zzfrld, zztms 
    9191      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zi1, zi2, zmasst, zpresh 
    92  
    9392      !!------------------------------------------------------------------- 
    94  
    95 !!bug 
    96 !!    u_oce(:,:) = 0.e0 
    97 !!    v_oce(:,:) = 0.e0 
    98 !!    write(*,*) 'rhg min, max u & v', maxval(u_oce), minval(u_oce), maxval(v_oce), minval(v_oce) 
    99 !!bug 
    10093       
    10194      !  Store initial velocities 
     
    107100      zztms(:,1:jpj) = tms(:,:)   ;    zzfrld(:,1:jpj) = frld(:,:) 
    108101      zu0(:,1:jpj) = u_ice(:,:)   ;    zv0(:,1:jpj) = v_ice(:,:) 
    109  
     102      ! 
    110103      zu_a(:,:)    = zu0(:,:)     ;   zv_a(:,:) = zv0(:,:) 
    111104      zu_n(:,:)    = zu0(:,:)     ;   zv_n(:,:) = zv0(:,:) 
     
    127120      zviszeta(:,:) = 0.e0 
    128121      zviseta (:,:) = 0.e0 
    129  
    130 !i    zviszeta(:,0    ) = 0.e0    ;    zviseta(:,0    ) = 0.e0 
    131 !i    zviszeta(:,jpj  ) = 0.e0    ;    zviseta(:,jpj  ) = 0.e0 
    132 !i    zviszeta(:,jpj+1) = 0.e0    ;    zviseta(:,jpj+1) = 0.e0 
    133122 
    134123 
     
    373362            END DO 
    374363         END DO 
    375  
    376          ! GAUSS-SEIDEL method 
     364         !  
    377365         !                                                      ! ================ ! 
    378 iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! 
     366iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! GAUSS-SEIDEL method 
    379367            !                                                   ! ================ ! 
    380368!CDIR NOVERRCHK 
  • branches/dev_1784_EVP/NEMO/LIM_SRC_2/limsbc_2.F90

    r2046 r2095  
    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 
    9    !!           09-05 (G.Garric) addition of the lim2_evp case 
     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 
    1010   !!---------------------------------------------------------------------- 
    1111#if defined key_lim2 
    1212   !!---------------------------------------------------------------------- 
    1313   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    14    !!---------------------------------------------------------------------- 
    1514   !!---------------------------------------------------------------------- 
    1615   !!   lim_sbc_2  : flux at the ice / ocean interface 
     
    1817   USE par_oce          ! ocean parameters 
    1918   USE dom_oce          ! ocean domain 
    20    USE sbc_ice          ! surface boundary condition 
    21    USE sbc_oce          ! surface boundary condition 
     19   USE sbc_ice          ! surface boundary condition: ice 
     20   USE sbc_oce          ! surface boundary condition: ocean 
    2221   USE phycst           ! physical constants 
    23    USE ice_2            ! LIM sea-ice variables 
    24  
    25    USE lbclnk           ! ocean lateral boundary condition 
     22   USE ice_2            ! LIM2: ice variables 
     23 
     24   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges 
    2625   USE in_out_manager   ! I/O manager 
    2726   USE diaar5, ONLY :   lk_diaar5 
    28    USE iom              !  
     27   USE iom              ! IOM library 
    2928   USE albedo           ! albedo parameters 
    3029   USE prtctl           ! Print control 
     
    3433   PRIVATE 
    3534 
    36    PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2 
    37  
    38    REAL(wp)  ::   epsi16 = 1.e-16  ! constant values 
    39    REAL(wp)  ::   rzero  = 0.e0     
    40    REAL(wp)  ::   rone   = 1.e0 
    41    REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r 
    42    REAL(wp), DIMENSION(jpi,jpj)  ::   sice_r 
     35   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2 
     36 
     37   REAL(wp)  ::   r1_rdtice = 1.e0 / rdt_ice   ! 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 
    4343 
    4444   !! * Substitutions 
    4545#  include "vectopt_loop_substitute.h90" 
    4646   !!---------------------------------------------------------------------- 
    47    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     47   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4848   !! $Id$ 
    4949   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5050   !!---------------------------------------------------------------------- 
    51  
    5251CONTAINS 
    5352 
     
    7978      !! 
    8079      INTEGER  ::   ji, jj           ! dummy loop indices 
    81       INTEGER  ::   ifvt, i1mfr, idfr               ! some switches 
    82       INTEGER  ::   iflt, ial, iadv, ifral, ifrdv 
    83       INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers 
    84       REAL(wp) ::   zrdtir           ! 1. / rdt_ice 
    85       REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux 
    86       REAL(wp) ::   zinda            ! switch for testing the values of ice concentration 
    87       REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface 
    88       REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface 
    89       REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points 
    90       REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points 
    91 !!!      REAL(wp) ::   zutaui , zvtaui  ! lead fraction at U- & V-points 
    92       REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity 
    93 ! interface 2D --> 3D 
    94       REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb     ! albedo of ice under overcast sky 
    95       REAL(wp), DIMENSION(jpi,jpj,1) ::   zalbp    ! albedo of ice under clear sky 
    96       REAL(wp) ::   zsang, zmod, zztmp, zfm 
    97       REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! component of ocean stress below sea-ice at I-point 
    98       REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi           ! module    of ocean stress below sea-ice at I-point 
    99       REAL(wp), DIMENSION(jpi,jpj) ::   zqnsoce          ! save qns before its modification by ice model 
    100  
     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 
    10190      !!--------------------------------------------------------------------- 
    10291      
    103       zrdtir = 1. / rdt_ice 
    10492       
    10593      IF( kt == nit000 ) THEN 
     
    10795         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition' 
    10896         IF(lwp) WRITE(numout,*) '~~~~~~~~~   ' 
    109  
     97         ! 
     98         r1_rdtice = 1. / rdt_ice 
     99         ! 
    110100         soce_r(:,:) = soce 
    111101         sice_r(:,:) = sice 
    112102         ! 
    113          IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN 
    114             !                                        ! ======================= 
    115             !                                        !  ORCA_R2 configuration 
    116             !                                        ! ======================= 
     103         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN     !  ORCA_R2 configuration 
    117104            ii0 = 145   ;   ii1 = 180        ! Baltic Sea 
    118105            ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 
     
    177164!!$ 
    178165 
    179             !   computation the solar flux at ocean surface 
     166            ! solar flux at ocean surface 
    180167#if defined key_coupled  
    181168            zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) ) 
     
    183170            zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj) 
    184171#endif             
    185             !  computation the non solar heat flux at ocean surface 
     172            ! non solar heat flux at ocean surface 
    186173            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads 
    187                &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            & 
    188                &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    & 
    189                &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir 
    190  
     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            ! 
    191178            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ??? 
    192              
     179            ! 
    193180            qsr  (ji,jj) = zqsr                                          ! solar heat flux  
    194181            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux 
     
    196183      END DO 
    197184 
    198       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )       
    199       CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )       
    200       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(:,:) ) ) 
    201188 
    202189      !------------------------------------------! 
    203190      !      mass flux at the ocean surface      ! 
    204191      !------------------------------------------! 
    205  
    206 !!gm 
    207 !!gm CAUTION    
    208 !!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm 
    209 !!gm 
    210192      DO jj = 1, jpj 
    211193         DO ji = 1, jpi 
    212              
    213194#if defined key_coupled 
    214           zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  
    215              &   + 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  
    216198#else 
    217 !!$            !  computing freshwater exchanges at the ice/ocean interface 
    218 !!$            zpme = - evap(ji,jj)    *   frld(ji,jj)           &   !  evaporation over oceanic fraction 
    219 !!$               &   + tprecip(ji,jj)                           &   !  total precipitation 
    220 !!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice 
    221 !!$               &   - rdmsnif(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting  
    222             !  computing freshwater exchanges at the ice/ocean interface 
    223             zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction  
    224                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean 
    225                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step 
    226                &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting  
    227                !                                                   !  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  
    228204#endif             
    229  
    230             !  computing salt exchanges at the ice/ocean interface 
     205            ! salt exchanges at the ice/ocean interface 
    231206            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir )  
    232              
    233             !  converting the salt flux from ice to a freshwater flux from ocean 
     207            ! 
     208            ! convert the salt flux from ice into a freshwater flux from ocean 
    234209            zfm  = zfons / ( sss_m(ji,jj) + epsi16 ) 
    235              
     210            ! 
    236211            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution) 
    237212            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution) 
    238  
    239213         END DO 
    240214      END DO 
    241  
     215      ! 
    242216      IF( lk_diaar5 ) THEN 
    243217         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * zrdtir ) 
     
    278252               ! ... components of ice-ocean stress at U and V-points  (from I-point values) 
    279253#if defined key_lim2_vp 
    280                zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) ) 
     254               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )      ! VP rheology 
    281255               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) ) 
    282256#else 
    283                zutau  = ztio_u(ji,jj) 
     257               zutau  = ztio_u(ji,jj)                                      ! EVP rheology 
    284258               zvtau  = ztio_v(ji,jj) 
    285259#endif 
     
    297271            END DO 
    298272         END DO 
    299  
    300          ! boundary condition on the stress (utau,vtau,taum) 
    301          CALL lbc_lnk( utau, 'U', -1. ) 
    302          CALL lbc_lnk( vtau, 'V', -1. ) 
     273         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition  
    303274         CALL lbc_lnk( taum, 'T',  1. ) 
    304275 
    305276      ENDIF 
    306277 
     278      IF( lk_cpl ) THEN            
    307279      !-----------------------------------------------! 
    308280      !   Coupling variables                          ! 
    309281      !-----------------------------------------------! 
    310  
    311       IF ( lk_cpl ) THEN            
    312          ! Ice surface temperature  
    313          tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature        
    314          ! Computation of snow/ice and ocean albedo 
     282         tn_ice(:,:,1) = sist(:,:)           ! sea-ice surface temperature        
     283         !                                   ! snow/ice and ocean albedo 
    315284         CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb ) 
    316285         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys) 
     
    325294         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ') 
    326295      ENDIF  
    327     
    328     END SUBROUTINE lim_sbc_2 
     296      ! 
     297   END SUBROUTINE lim_sbc_2 
    329298 
    330299#else 
  • branches/dev_1784_EVP/NEMO/LIM_SRC_2/limtrp_2.F90

    r2046 r2095  
    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 
    69       !!   3.3  !  09-05  (G.Garric) addition of the lim2_evp case 
    7062      !!--------------------------------------------------------------------- 
    7163      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    72  
    73       INTEGER  ::   ji, jj, jk,   &  ! dummy loop indices 
    74          &          initad           ! number of sub-timestep for the advection 
    75  
    76       REAL(wp) ::  &                               
    77          zindb  ,  & 
    78          zacrith, & 
    79          zindsn , & 
    80          zindic , & 
    81          zusvosn, & 
    82          zusvoic, & 
    83          zignm  , & 
    84          zindhe , & 
    85          zvbord , & 
    86          zcfl   , & 
    87          zusnit , & 
    88          zrtt, ztsn, ztic1, ztic2 
    89  
    90       REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace 
    91          zui_u , zvi_v , zsm   ,         & 
    92          zs0ice, zs0sn , zs0a  ,         & 
    93          zs0c0 , zs0c1 , zs0c2 ,         & 
    94          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   !  -      - 
    9575      !--------------------------------------------------------------------- 
    9676 
     
    10989         zvbord = 1.0 + ( 1.0 - bound ) 
    11090#if defined key_lim2_vp 
    111          DO jj = 1, jpjm1 
     91         DO jj = 1, jpjm1     ! VP rheology: ice (u,v) at I-point 
    11292            DO ji = 1, jpim1   ! NO vector opt. 
    113                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 ) ) 
    114                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 ) ) 
    11595            END DO 
    11696         END DO 
    117          ! Lateral boundary conditions on zui_u, zvi_v 
    118          CALL lbc_lnk( zui_u, 'U', -1. ) 
    119          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  
    12098#else 
    121         zui_u(:,:)=u_ice(:,:) 
    122         zvi_v(:,:)=v_ice(:,:) 
     99        zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points 
     100        zvi_v(:,:) = v_ice(:,:) 
    123101#endif 
    124102 
     
    129107         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) ) 
    130108 
    131          IF (lk_mpp ) CALL mpp_max(zcfl) 
    132  
    133          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 
    134112 
    135113         ! content of properties 
    136114         ! --------------------- 
    137          zs0sn (:,:) =  hsnm(:,:) * area(:,:)                 ! Snow volume. 
    138          zs0ice(:,:) =  hicm (:,:) * area(:,:)                ! Ice volume. 
    139          zs0a  (:,:) =  ( 1.0 - frld(:,:) ) * area(:,:)       ! Surface covered by ice. 
    140          zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn(:,:)   ! Heat content of the snow layer. 
    141          zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer. 
    142          zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer. 
    143          zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a(:,:)    ! Heat reservoir for brine pockets. 
    144           
     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. 
    145122  
    146123         ! Advection  
     
    150127         zusnit = 1.0 / REAL( initad )  
    151128          
    152          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 
    153132            DO jk = 1,initad 
    154133               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) 
     
    237216         ! Up-dating and limitation of sea ice properties after transport. 
    238217         DO jj = 1, jpj 
    239 !!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq! 
    240218            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH 
    241219            DO ji = 1, jpi 
     
    278256            END DO 
    279257         END DO 
    280           
     258         ! 
    281259      ENDIF 
    282        
     260      ! 
    283261   END SUBROUTINE lim_trp_2 
    284262 
     
    294272      !! 
    295273      !! ** input   :   Namelist namicetrp 
    296       !! 
    297       !! history : 
    298       !!   2.0  !  03-08 (C. Ethe)  Original code 
    299274      !!------------------------------------------------------------------- 
    300275      NAMELIST/namicetrp/ bound 
    301276      !!------------------------------------------------------------------- 
    302  
     277      ! 
    303278      ! Read Namelist namicetrp 
    304279      REWIND ( numnam_ice ) 
     
    310285         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound 
    311286      ENDIF 
    312              
     287      !      
    313288   END SUBROUTINE lim_trp_init_2 
    314289 
  • branches/dev_1784_EVP/NEMO/LIM_SRC_3/limrhg.F90

    r2046 r2095  
    44   !!   Ice rheology : sea ice rheology 
    55   !!====================================================================== 
    6    !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code 
    7    !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3 
     6   !! History :  LIM  !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code 
     7   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM 3 
    88   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy  
    9    !!             -   !  2009-05  (G.Garric) addition of the lim2_evp cas 
     9   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas 
    1010   !!---------------------------------------------------------------------- 
    1111#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp ) 
     
    1515   !!   lim_rhg   : computes ice velocities 
    1616   !!---------------------------------------------------------------------- 
    17    !! * Modules used 
    18    USE phycst 
    19    USE par_oce 
    20    USE dom_oce 
    21    USE sbc_oce         ! Surface boundary condition: ocean fields 
    22    USE sbc_ice         ! Surface boundary condition: ice fields 
    23    USE lbclnk 
    24    USE lib_mpp 
    25    USE in_out_manager  ! I/O manager 
    26    USE limitd_me 
    27    USE prtctl          ! Print control 
     17   USE phycst           ! physical constants 
     18   USE par_oce          ! ocean parameters 
     19   USE dom_oce          ! ocean domain 
     20   USE sbc_oce          ! Surface boundary condition: ocean fields 
     21   USE sbc_ice          ! Surface boundary condition: ice fields 
     22   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     23   USE lib_mpp          ! MPP library 
     24   USE in_out_manager   ! I/O manager 
     25   USE limitd_me        ! LIM3:  
     26   USE prtctl           ! control print 
    2827#if defined key_lim3 
    29    USE dom_ice 
    30    USE ice 
    31    USE iceini 
     28   USE ice              ! LIM3: ice variables 
     29   USE dom_ice          ! LIM3: ice domain 
     30   USE iceini           ! LIM3: ice initialisation 
    3231#endif 
    3332#if defined key_lim2 && ! defined key_lim2_vp 
    34    USE dom_ice_2 
    35    USE ice_2 
    36    USE iceini_2 
     33   USE ice_2            ! LIM2: ice variables 
     34   USE dom_ice_2        ! LIM2: ice domain 
     35   USE iceini_2         ! LIM2: ice initialisation 
    3736#endif 
    3837 
     
    4039   PRIVATE 
    4140 
    42    !! * Routine accessibility 
    43    PUBLIC lim_rhg  ! routine called by lim_dyn 
    44  
    45    !! * Substitutions 
    46 #  include "vectopt_loop_substitute.h90" 
     41   PUBLIC   lim_rhg   ! routine called by lim_dyn module 
    4742 
    4843   !! * Module variables 
     
    5045      rzero   = 0.e0   ,  & 
    5146      rone    = 1.e0 
     47 
     48   !! * Substitutions 
     49#  include "vectopt_loop_substitute.h90" 
    5250   !!---------------------------------------------------------------------- 
    53    !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)  
     51   !! NEMO/LIM-3 3.3,  UCL-LOCEAN-IPSL (2010)  
    5452   !! $Id$ 
    5553   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5654   !!---------------------------------------------------------------------- 
    57  
    5855CONTAINS 
    5956 
    6057   SUBROUTINE lim_rhg( k_j1, k_jpj ) 
    61  
    6258      !!------------------------------------------------------------------- 
    6359      !!                 ***  SUBROUTINE lim_rhg  *** 
     
    107103      !! 
    108104      !! ** References : Hunke and Dukowicz, JPO97 
    109       !!                 Bouillon et al., 08, in prep (update this when 
    110       !!                 published) 
    111       !!                 Vancoppenolle et al., OM08 
     105      !!                 Bouillon et al., 2009, Ocean. Modelling, 27, 174-184. 
     106      !!                 Vancoppenolle et al. 2009, Ocean Modelling, 27, 33-53. 
     107      !!------------------------------------------------------------------- 
     108      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation 
     109      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation 
    112110      !! 
    113       !!------------------------------------------------------------------- 
    114       ! * Arguments 
    115       ! 
    116       INTEGER, INTENT(in) :: & 
    117          k_j1 ,                      & !: southern j-index for ice computation 
    118          k_jpj                         !: northern j-index for ice computation 
    119  
    120       ! * Local variables 
    121       INTEGER ::   ji, jj              !: dummy loop indices 
    122  
    123       INTEGER  :: & 
    124          jter                          !: temporary integers 
    125  
    126       CHARACTER (len=50) ::   charout 
    127  
    128       REAL(wp) :: & 
    129          zt11, zt12, zt21, zt22,     & !: temporary scalars 
    130          ztagnx, ztagny,             & !: wind stress on U/V points                        
    131          delta                         ! 
    132  
    133       REAL(wp) :: & 
    134          za,                         & !: 
    135          zstms,                      & !: temporary scalar for ice strength 
    136          zsang,                      & !: temporary scalar for coriolis term 
    137          zmask                         !: mask for the computation of ice mass 
    138  
    139       REAL(wp),DIMENSION(jpi,jpj) :: & 
    140          zpresh        ,             & !: temporary array for ice strength 
    141          zpreshc       ,             & !: Ice strength on grid cell corners (zpreshc) 
    142          zfrld1, zfrld2,             & !: lead fraction on U/V points                                     
    143          zmass1, zmass2,             & !: ice/snow mass on U/V points                                     
    144          zcorl1, zcorl2,             & !: coriolis parameter on U/V points 
    145          za1ct, za2ct  ,             & !: temporary arrays 
    146          zc1           ,             & !: ice mass 
    147          zusw          ,             & !: temporary weight for the computation 
    148                                 !: of ice strength 
    149          u_oce1, v_oce1,             & !: ocean u/v component on U points                            
    150          u_oce2, v_oce2,             & !: ocean u/v component on V points 
    151          u_ice2,                     & !: ice u component on V point 
    152          v_ice1                        !: ice v component on U point 
     111      INTEGER  ::   ji, jj   ! dummy loop indices 
     112      INTEGER  ::   jter     ! local integers 
     113      CHARACTER (len=50) ::   charout   ! local character 
     114      REAL(wp) ::   zt11, zt12, zt21, zt22   ! local scalars 
     115      REAL(wp) ::   ztagnx, ztagny, delta    !   -      - 
     116      REAL(wp) ::   za, zstms, zsang, zmask  !   -      - 
     117      REAL(wp) ::   zresm, zindb, zdummy     !   -    - 
     118 
     119      REAL(wp),DIMENSION(jpi,jpj) ::   zpresh , zfrld1, zmass1, zcorl1, za1ct    ! 2D workspace 
     120      REAL(wp),DIMENSION(jpi,jpj) ::   zpreshc, zfrld2, zmass2, zcorl2, za2ct    !  -      - 
     121      REAL(wp),DIMENSION(jpi,jpj) ::   u_oce1, v_oce1, u_ice2, zc1               !  -      -                           
     122      REAL(wp),DIMENSION(jpi,jpj) ::   u_oce2, v_oce2, v_ice1, zusw              !  -      - 
     123      REAL(wp),DIMENSION(jpi,jpj) ::   zf1, zf2                                  !  -      - 
    153124 
    154125      REAL(wp) :: & 
     
    166137         sigma1, sigma2                !: internal ice stress 
    167138 
    168       REAL(wp),DIMENSION(jpi,jpj) :: & 
    169          zf1, zf2                      !: arrays for internal stresses 
    170139 
    171140      REAL(wp),DIMENSION(jpi,jpj) :: & 
     
    177146         zs12                          !: Non-diagonal stress tensor component zs12 
    178147 
    179       REAL(wp) :: & 
    180          zresm            ,          & !: Maximal error on ice velocity 
    181          zindb            ,          & !: ice (1) or not (0)       
    182          zdummy                        !: dummy argument 
    183  
    184       REAL(wp),DIMENSION(jpi,jpj) :: & 
    185          zu_ice           ,          & !: Ice velocity on previous time step 
    186          zv_ice           ,          & 
    187          zresr                         !: Local error on velocity 
     148 
     149      REAL(wp),DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr   ! 
     150      !!------------------------------------------------------------------- 
    188151 
    189152#if  defined key_lim2 && ! defined key_lim2_vp 
Note: See TracChangeset for help on using the changeset viewer.