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

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/LIM_SRC_2
Files:
8 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 
Note: See TracChangeset for help on using the changeset viewer.