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 888 for trunk/NEMO/LIM_SRC_2 – NEMO

Changeset 888 for trunk/NEMO/LIM_SRC_2


Ignore:
Timestamp:
2008-04-11T19:05:03+02:00 (16 years ago)
Author:
ctlod
Message:

merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113

Location:
trunk/NEMO/LIM_SRC_2
Files:
1 added
2 deleted
21 edited

Legend:

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

    r823 r888  
    11MODULE dom_ice_2 
    2 #if defined key_lim2 
    32   !!====================================================================== 
    43   !!                   ***  MODULE  dom_ice  *** 
     
    76   !! History :   2.0  !  03-08  (C. Ethe)  Free form and module 
    87   !!---------------------------------------------------------------------- 
    9  
     8#if defined key_lim2 
    109   !!---------------------------------------------------------------------- 
    1110   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    12    !! $Header$ 
     11   !! $ Id: $ 
    1312   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    1413   !!---------------------------------------------------------------------- 
  • trunk/NEMO/LIM_SRC_2/ice_2.F90

    r823 r888  
    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   !!---------------------------------------------------------------------- 
    68#if defined key_lim2 
    79   !!---------------------------------------------------------------------- 
    810   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    911   !!---------------------------------------------------------------------- 
    10    !! History : 
    11    !!   2.0  !  03-08  (C. Ethe)  F90: Free form and module 
    12    !!---------------------------------------------------------------------- 
    13    !!  LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    14    !! $Header$ 
    15    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
     12   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006) 
     13   !! $ Id: $ 
     14   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    1615   !!---------------------------------------------------------------------- 
    1716   !! * Modules used 
     
    2120   PRIVATE 
    2221 
    23    !! * Share Module variables 
    24    INTEGER , PUBLIC ::   & !!: ** ice-dynamic namelist (namicedyn) ** 
    25       nbiter = 1      ,  &  !: number of sub-time steps for relaxation 
    26       nbitdr = 250          !: maximum number of iterations for relaxation 
     22   !!* ice-dynamic namelist (namicedyn) * 
     23   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation 
     24   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation 
     25   REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic 
     26   REAL(wp), PUBLIC ::   alpha  = 0.5       !: coefficient for semi-implicit coriolis 
     27   REAL(wp), PUBLIC ::   dm     = 0.6e+03   !: diffusion constant for dynamics 
     28   REAL(wp), PUBLIC ::   om     = 0.5       !: relaxation constant 
     29   REAL(wp), PUBLIC ::   resl   = 5.0e-05   !: maximum value for the residual of relaxation 
     30   REAL(wp), PUBLIC ::   cw     = 5.0e-03   !: drag coefficient for oceanic stress 
     31   REAL(wp), PUBLIC ::   angvg  = 0.e0      !: turning angle for oceanic stress 
     32   REAL(wp), PUBLIC ::   pstar  = 1.0e+04   !: first bulk-rheology parameter 
     33   REAL(wp), PUBLIC ::   c_rhg  = 20.e0     !: second bulk-rhelogy parameter 
     34   REAL(wp), PUBLIC ::   etamn  = 0.e+07    !: minimun value for viscosity 
     35   REAL(wp), PUBLIC ::   creepl = 2.e-08    !: creep limit 
     36   REAL(wp), PUBLIC ::   ecc    = 2.e0      !: eccentricity of the elliptical yield curve 
     37   REAL(wp), PUBLIC ::   ahi0   = 350.e0    !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    2738 
    28    REAL(wp), PUBLIC ::   & !!: ** ice-dynamic namelist (namicedyn) ** 
    29       epsd   = 1.0e-20,  &  !: tolerance parameter for dynamic 
    30       alpha  = 0.5    ,  &  !: coefficient for semi-implicit coriolis 
    31       dm     = 0.6e+03,  &  !: diffusion constant for dynamics 
    32       om     = 0.5    ,  &  !: relaxation constant 
    33       resl   = 5.0e-05,  &  !: maximum value for the residual of relaxation 
    34       cw     = 5.0e-03,  &  !: drag coefficient for oceanic stress 
    35       angvg  = 0.e0   ,  &  !: turning angle for oceanic stress 
    36       pstar  = 1.0e+04,  &  !: first bulk-rheology parameter 
    37       c_rhg  = 20.e0  ,  &  !: second bulk-rhelogy parameter 
    38       etamn  = 0.e+07,   &  !: minimun value for viscosity 
    39       creepl = 2.e-08,   &  !: creep limit 
    40       ecc    = 2.e0   ,  &  !: eccentricity of the elliptical yield curve 
    41       ahi0   = 350.e0       !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
     39   REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc ) 
     40   REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw 
     41   REAL(wp), PUBLIC ::   sangvg, cangvg     !: sin and cos of the turning angle for ocean stress 
     42   REAL(wp), PUBLIC ::   pstarh             !: pstar / 2.0 
    4243 
    43    REAL(wp), PUBLIC ::   &  !: 
    44       usecc2          ,  &  !:  = 1.0 / ( ecc * ecc ) 
    45       rhoco           ,  &  !: = rau0 * cw 
    46       sangvg, cangvg  ,  &  !: sin and cos of the turning angle for ocean stress 
    47       pstarh                !: pstar / 2.0 
     44   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
     45   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
     46   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses 
     47   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity 
    4848 
    49    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &  !: 
    50       u_oce, v_oce,      &  !: surface ocean velocity used in ice dynamics 
    51       ahiu , ahiv ,      &  !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
    52       pahu , pahv ,      &  !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
    53       hsnm , hicm ,      &  !: mean snow and ice thicknesses 
    54       ust2s                 !: friction velocity 
     49   !!* diagnostic quantities 
     50!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   firic         !: IR flux over the ice (only used for outputs) 
     51!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fcsic         !: Sensible heat flux over the ice (only used for outputs) 
     52!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fleic         !: Latent heat flux over the ice (only used for outputs) 
     53!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qlatic        !: latent flux (only used for outputs) 
     54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvosif       !: Variation of volume at surface (only used for outputs) 
     55   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvobif       !: Variation of ice volume at the bottom ice (only used for outputs) 
     56   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdvolif       !: Total variation of ice volume (only used for outputs) 
     57   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvonif       !: Lateral Variation of ice volume (only used for outputs) 
    5558 
    56    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &  !: 
    57         sst_ini,         &  !: sst read from a file for ice model initialization  
    58         sss_ini             !: sss read from a file for ice model initialization  
     59   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin ??? degree ??? I don't know) 
     60   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS 
     61   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness 
     62   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnif         !: Snow thickness 
     63   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicifp        !: Ice production/melting 
     64   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   frld          !: Leads fraction = 1-a/totalarea 
     65   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   phicif        !: ice thickness  at previous time  
     66   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pfrld         !: Leads fraction at previous time   
     67   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qstoif        !: Energy stored in the brine pockets 
     68   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fbif          !: Heat flux at the ice base 
     69   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmsnif       !: Variation of snow mass 
     70   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass 
     71   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean) 
     72   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing  
     73   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean 
     74   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2) 
     75   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget 
     76   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice 
     77   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max heat contained in brine pockets (?) 
     78   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?) 
     79   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?) 
     80   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy in case of toral lateral ablation (?) 
     81   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice 
    5982 
    60    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    61       firic  ,   &  !: IR flux over the ice (only used for outputs) 
    62       fcsic  ,   &  !: Sensible heat flux over the ice (only used for outputs) 
    63       fleic  ,   &  !: Latent heat flux over the ice (only used for outputs) 
    64       qlatic ,   &  !: latent flux 
    65       rdvosif,   &  !: Variation of volume at surface (only used for outputs) 
    66       rdvobif,   &  !: Variation of ice volume at the bottom ice (only used for outputs) 
    67       fdvolif,   &  !: Total variation of ice volume (only used for outputs) 
    68       rdvonif,   &  !: Lateral Variation of ice volume (only used for outputs) 
    69       sist   ,   &  !: Sea-Ice Surface Temperature (Kelvin ??? degree ??? I don't know) 
    70       tfu    ,   &  !: Melting point temperature of sea water 
    71       hsnif  ,   &  !: Snow thickness 
    72       hicif  ,   &  !: Ice thickness 
    73       hicifp ,   &  !: Ice production/melting 
    74       frld   ,   &  !: Leads fraction = 1-a/totalarea 
    75       phicif ,   &  !: ice thickness  at previous time  
    76       pfrld  ,   &  !: Leads fraction at previous time   
    77       qstoif ,   &  !: Energy stored in the brine pockets 
    78       fbif   ,   &  !: Heat flux at the ice base 
    79       rdmsnif,   &  !: Variation of snow mass 
    80       rdmicif,   &  !: Variation of ice mass 
    81       qldif  ,   &  !: heat balance of the lead (or of the open ocean) 
    82       qcmif  ,   &  !: Energy needed to bring the ocean surface layer until its freezing  
    83       fdtcn  ,   &  !: net downward heat flux from the ice to the ocean 
    84       qdtcn  ,   &  !: energy from the ice to the ocean 
    85       !             !  point (at a factor 2) 
    86       thcm   ,   &  !: part of the solar energy used in the lead heat budget 
    87       fstric ,   &  !: Solar flux transmitted trough the ice 
    88       ffltbif,   &  !: Array linked with the max heat contained in brine pockets (?) 
    89       fscmbq ,   &  !: Linked with the solar flux below the ice (?) 
    90       fsbbq  ,   &  !: Also linked with the solar flux below the ice (?) 
    91       qfvbq  ,   &  !: Array used to store energy in case of toral lateral ablation (?) 
    92       dmgwi         !: Variation of the mass of snow ice 
     83   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   albege        !: Albedo of the snow or ice (only for outputs) 
     84   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   albecn        !: Albedo of the ocean (only for outputs) 
     85   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tauc          !: Cloud optical depth 
    9386 
    94    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    95       albege ,   &  !: Albedo of the snow or ice (only for outputs) 
    96       albecn ,   &  !: Albedo of the ocean (only for outputs) 
    97       tauc   ,   &  !: Cloud optical depth 
    98       sdvt          !: u*^2/(Stress/density) 
     87   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ui_ice, vi_ice   !: two components of the ice   velocity at I-point (m/s) 
     88   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ui_oce, vi_oce   !: two components of the ocean velocity at I-point (m/s) 
    9989 
     90   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpsmax)     ::   scal0   !: ??? 
     91   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer 
    10092 
    101    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    102       u_ice, v_ice,   &  !: two components of the ice velocity (m/s) 
    103       tio_u, tio_v       !: two components of the ice-ocean stress (N/m2) 
     93!! REAL(wp), DIMENSION(jpi,jpj,0:jpkmax+1) ::   reslum        !: Relative absorption of solar radiation in each ocean level 
    10494 
    105    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpsmax) ::   &  !: 
    106       scal0              !: ??? 
    107  
    108    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   &  !: 
    109       tbif          !: Temperature inside the ice/snow layer 
    110  
    111    REAL(wp), DIMENSION(jpi,jpj,0:jpkmax+1) ::    &  !: 
    112       reslum        !: Relative absorption of solar radiation in each ocean level 
    113  
    114    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    115          sxice, syice, sxxice, syyice, sxyice,      &  !: moments for advection 
    116          sxsn,  sysn,  sxxsn,  syysn,  sxysn,       &  !: 
    117          sxa,   sya,   sxxa,   syya,   sxya,        &  !: 
    118          sxc0,  syc0,  sxxc0,  syyc0,  sxyc0,       &  !: 
    119          sxc1,  syc1,  sxxc1,  syyc1,  sxyc1,       &  !: 
    120          sxc2,  syc2,  sxxc2,  syyc2,  sxyc2,       &  !: 
    121          sxst,  syst,  sxxst,  syyst,  sxyst           !: 
     95   !!* moment used in the advection scheme 
     96   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume 
     97   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume 
     98   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxa,   sya,   sxxa,   syya,   sxya     !: for ice cover area 
     99   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc0,  syc0,  sxxc0,  syyc0,  sxyc0    !: for heat content of snow 
     100   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc1,  syc1,  sxxc1,  syyc1,  sxyc1    !: for heat content of 1st ice layer 
     101   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc2,  syc2,  sxxc2,  syyc2,  sxyc2    !: for heat content of 2nd ice layer 
     102   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxst,  syst,  sxxst,  syyst,  sxyst    !: for heat content of brine pockets 
    122103 
    123104#else 
  • trunk/NEMO/LIM_SRC_2/iceini_2.F90

    r823 r888  
    1717   USE dom_oce 
    1818   USE dom_ice_2 
    19    USE in_out_manager 
    2019   USE ice_oce         ! ice variables 
    21    USE flx_oce 
     20   USE sbc_oce         ! surface boundary condition: ocean 
     21   USE sbc_ice         ! surface boundary condition: ice 
    2222   USE phycst          ! Define parameters for the routines 
    2323   USE ocfzpt 
     
    2727   USE limrst_2    
    2828   USE ini1d           ! initialization of the 1D configuration 
     29   USE in_out_manager 
    2930       
    3031   IMPLICIT NONE 
     
    4041   !!---------------------------------------------------------------------- 
    4142   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    42    !! $Header$  
     43   !! $ Id: $  
    4344   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    4445   !!---------------------------------------------------------------------- 
     
    6263                  
    6364      ! Louvain la Neuve Ice model 
    64       IF( nacc == 1 ) THEN 
    65           dtsd2   = nfice * rdtmin * 0.5 
    66           rdt_ice = nfice * rdtmin 
    67       ELSE 
    68           dtsd2   = nfice * rdt * 0.5 
    69           rdt_ice = nfice * rdt 
    70       ENDIF 
     65      dtsd2   = nn_fsbc * rdttra(1) * 0.5 
     66      rdt_ice = nn_fsbc * rdttra(1) 
    7167 
    7268      CALL lim_msh_2                  ! ice mesh initialization 
  • trunk/NEMO/LIM_SRC_2/limadv_2.F90

    r823 r888  
    3333   !!---------------------------------------------------------------------- 
    3434   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    35    !! $Header$  
     35   !! $ Id: $  
    3636   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3737   !!---------------------------------------------------------------------- 
  • trunk/NEMO/LIM_SRC_2/limdia_2.F90

    r823 r888  
    1919   USE par_ice_2       ! ice parameters 
    2020   USE ice_oce         ! ice variables 
     21   USE sbc_oce         ! surface boundary condition variables 
    2122   USE daymod          ! 
    2223   USE dom_ice_2       ! 
     
    2829   PRIVATE 
    2930 
    30    PUBLIC               lim_dia_2          ! called by ice_step 
     31   PUBLIC               lim_dia_2          ! called by sbc_ice_lim_2 
    3132   INTEGER, PUBLIC ::   ntmoy   = 1 ,   &  !: instantaneous values of ice evolution or averaging ntmoy 
    3233      &                 ninfo   = 1        !: frequency of ouputs on file ice_evolu in case of averaging 
     
    5859   !!---------------------------------------------------------------------- 
    5960   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    60    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdia.F90,v 1.9 2007/06/29 17:03:12 opalod Exp $  
     61   !! $ Id: $ 
    6162   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    6263   !!---------------------------------------------------------------------- 
     
    8788 
    8889      nv = 1  
    89       vinfor(nv) = REAL( kt + nfice - 1 ) 
     90      vinfor(nv) = REAL( kt + nn_fsbc - 1 ) 
    9091      nv = nv + 1 
    9192      vinfor(nv) = nyear 
     
    107108               zicevol = zarea   * hicif(ji,jj) 
    108109               zsnwvol = zarea   * hsnif(ji,jj) 
    109                zicespd = zicevol * ( u_ice(ji,jj) * u_ice(ji,jj)   & 
    110                   &                + v_ice(ji,jj) * v_ice(ji,jj) ) 
     110               zicespd = zicevol * ( ui_ice(ji,jj) * ui_ice(ji,jj)   & 
     111                  &                + vi_ice(ji,jj) * vi_ice(ji,jj) ) 
    111112               vinfor(nv+ 1) = vinfor(nv+ 1) + zarea 
    112113               vinfor(nv+ 3) = vinfor(nv+ 3) + zextent15 
     
    133134                zicevol = zarea   * hicif(ji,jj) 
    134135                zsnwvol = zarea   * hsnif(ji,jj) 
    135                 zicespd = zicevol * ( u_ice(ji,jj) * u_ice(ji,jj)   & 
    136                    &                + v_ice(ji,jj) * v_ice(ji,jj) ) 
     136                zicespd = zicevol * ( ui_ice(ji,jj) * ui_ice(ji,jj)   & 
     137                   &                + vi_ice(ji,jj) * vi_ice(ji,jj) ) 
    137138                vinfor(nv+ 1) = vinfor(nv+ 1) + zarea 
    138139                vinfor(nv+ 3) = vinfor(nv+ 3) + zextent15 
     
    154155     
    155156       ! oututs on file ice_evolu     
    156        IF( MOD( kt + nfice - 1, ninfo ) == 0 ) THEN 
     157       IF( MOD( kt + nn_fsbc - 1, ninfo ) == 0 ) THEN 
    157158          WRITE(numevo_ice,fmtw) ( titvar(jv), vinfom(jv)/naveg, jv = 1, nvinfo ) 
    158159          naveg = 0 
     
    227228 
    228229       ! Definition et Ecriture de l'entete : nombre d'enregistrements  
    229        ndeb   = ( nit000 - 1 + nfice - 1 ) / ninfo 
    230        IF( nit000 - 1 + nfice == 1 ) ndeb = -1 
    231  
    232        nferme = ( nitend + nfice - 1 ) / ninfo ! nit000 - 1 + nfice - 1 + nitend - nit000 + 1 
     230       ndeb   = ( nit000 - 1 + nn_fsbc - 1 ) / ninfo 
     231       IF( nit000 - 1 + nn_fsbc == 1 ) ndeb = -1 
     232 
     233       nferme = ( nitend + nn_fsbc - 1 ) / ninfo ! nit000 - 1 + nn_fsbc - 1 + nitend - nit000 + 1 
    233234       ntot   = nferme - ndeb 
    234235       ndeb   = ninfo * ( 1 + ndeb ) 
  • trunk/NEMO/LIM_SRC_2/limdmp_2.F90

    r823 r888  
    4040   !!---------------------------------------------------------------------- 
    4141   !!   LIM 2.0 , UCL-LOCEAN-IPSL  (2006) 
    42    !! $Header$ 
     42   !! $ Id: $ 
    4343   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4444   !!---------------------------------------------------------------------- 
  • trunk/NEMO/LIM_SRC_2/limdyn_2.F90

    r823 r888  
    44   !!   Sea-Ice dynamics :   
    55   !!====================================================================== 
     6   !! History :   1.0  !  01-04  (LIM)  Original code 
     7   !!             2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp 
     8   !!             2.0  !  03-08  (C. Ethe) add lim_dyn_init 
     9   !!             2.0  !  06-07  (G. Madec)  Surface module 
     10   !!--------------------------------------------------------------------- 
    611#if defined key_lim2 
    712   !!---------------------------------------------------------------------- 
     
    1116   !!    lim_dyn_init_2 : initialization and namelist read 
    1217   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    14    USE phycst 
    15    USE in_out_manager  ! I/O manager 
    16    USE dom_ice_2 
    17    USE dom_oce         ! ocean space and time domain 
    18    USE ice_2 
    19    USE ice_oce 
    20    USE iceini_2 
    21    USE limistate_2 
    22    USE limrhg_2        ! ice rheology 
    23    USE lbclnk 
    24    USE lib_mpp 
    25    USE prtctl          ! Print control 
     18   USE dom_oce        ! ocean space and time domain 
     19   USE sbc_oce        ! 
     20   USE phycst         ! 
     21   USE ice_2          ! 
     22   USE ice_oce        ! 
     23   USE dom_ice_2      ! 
     24   USE iceini_2       ! 
     25   USE limistate_2    ! 
     26   USE limrhg_2       ! ice rheology 
     27 
     28   USE lbclnk         ! 
     29   USE lib_mpp        ! 
     30   USE in_out_manager ! I/O manager 
     31   USE prtctl         ! Print control 
    2632 
    2733   IMPLICIT NONE 
    2834   PRIVATE 
    2935 
    30    !! * Accessibility 
    31    PUBLIC lim_dyn_2  ! routine called by ice_step 
     36   PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim 
    3237 
    3338   !! * Module variables 
    3439   REAL(wp)  ::  rone    = 1.e0   ! constant value 
    3540 
    36    !!---------------------------------------------------------------------- 
    37    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    38    !! $Header$  
    39    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     41#  include "vectopt_loop_substitute.h90" 
     42   !!---------------------------------------------------------------------- 
     43   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     44   !! $ Id: $ 
     45   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4046   !!---------------------------------------------------------------------- 
    4147 
     
    4652      !!               ***  ROUTINE lim_dyn_2  *** 
    4753      !!                
    48       !! ** Purpose :   compute ice velocity and ocean-ice stress 
     54      !! ** Purpose :   compute ice velocity and ocean-ice friction velocity 
    4955      !!                 
    5056      !! ** Method  :  
     
    5258      !! ** Action  : - Initialisation 
    5359      !!              - Call of the dynamic routine for each hemisphere 
    54       !!              - computation of the stress at the ocean surface          
     60      !!              - computation of the friction velocity at the sea-ice base 
    5561      !!              - treatment of the case if no ice dynamic 
    56       !! History : 
    57       !!   1.0  !  01-04  (LIM)  Original code 
    58       !!   2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp 
    5962      !!--------------------------------------------------------------------- 
    6063      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    61  
    62       INTEGER ::   ji, jj             ! dummy loop indices 
    63       INTEGER ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology 
    64       REAL(wp) ::   & 
    65          ztairx, ztairy,           &  ! tempory scalars 
    66          zsang , zmod,             & 
    67          ztglx , ztgly ,           & 
    68          zt11, zt12, zt21, zt22 ,  & 
    69          zustm, zsfrld, zsfrldm4,  & 
    70          zu_ice, zv_ice, ztair2 
    71       REAL(wp),DIMENSION(jpj) ::   & 
    72          zind,                     &  ! i-averaged indicator of sea-ice 
    73          zmsk                         ! i-averaged of tmask 
     64      !! 
     65      INTEGER  ::   ji, jj             ! dummy loop indices 
     66      INTEGER  ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology 
     67      REAL(wp) ::   zcoef              ! temporary scalar 
     68      REAL(wp), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice 
     69      REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask 
     70      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io   ! ice-ocean velocity 
    7471      !!--------------------------------------------------------------------- 
    7572 
    76       IF( kt == nit000  )   CALL lim_dyn_init_2   ! Initialization (first time-step only) 
     73      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only) 
    7774       
    78       IF ( ln_limdyn ) THEN 
    79  
     75      IF( ln_limdyn ) THEN 
     76         ! 
    8077         ! Mean ice and snow thicknesses.           
    8178         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:) 
    8279         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:) 
    83  
    84          u_oce(:,:)  = u_io(:,:) * tmu(:,:) 
    85          v_oce(:,:)  = v_io(:,:) * tmu(:,:) 
    86         
    87          !                                         ! Rheology (ice dynamics) 
    88          !                                         ! ======== 
     80         ! 
     81         !                                     ! Rheology (ice dynamics) 
     82         !                                     ! ======== 
    8983          
    9084         !  Define the j-limits where ice rheology is computed 
     
    9488            i_j1 = 1    
    9589            i_jpj = jpj 
    96             IF(ln_ctl)    THEN 
    97                CALL prt_ctl_info('lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj) 
    98             ENDIF 
     90            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    9991            CALL lim_rhg_2( i_j1, i_jpj ) 
    100  
     92            ! 
    10193         ELSE                                 ! optimization of the computational area 
    102  
     94            ! 
    10395            DO jj = 1, jpj 
    10496               zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line 
    10597               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line 
    106    !!i         write(numout,*) narea, 'limdyn' , jj, zind(jj), zmsk(jj) 
    107             END DO 
    108  
     98            END DO 
     99            ! 
    109100            IF( l_jeq ) THEN                     ! local domain include both hemisphere 
    110101               !                                 ! Rheology is computed in each hemisphere 
     
    118109               i_j1 = MAX( 1, i_j1-1 ) 
    119110               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    120      
     111               !  
    121112               CALL lim_rhg_2( i_j1, i_jpj ) 
    122      
     113               !  
    123114               ! Southern hemisphere 
    124115               i_j1  =  1  
     
    129120               i_jpj = MIN( jpj, i_jpj+2 ) 
    130121               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    131      
     122               !  
    132123               CALL lim_rhg_2( i_j1, i_jpj ) 
    133      
     124               !  
    134125            ELSE                                 ! local domain extends over one hemisphere only 
    135126               !                                 ! Rheology is computed only over the ice cover 
     
    148139     
    149140               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    150      
     141               !  
    151142               CALL lim_rhg_2( i_j1, i_jpj ) 
    152  
     143               ! 
    153144            ENDIF 
    154  
     145            ! 
    155146         ENDIF 
    156147 
    157          IF(ln_ctl)   THEN  
    158             CALL prt_ctl(tab2d_1=u_oce , clinfo1=' lim_dyn  : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce :') 
    159             CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 
    160          ENDIF 
    161           
    162          !                                         ! Ice-Ocean stress 
    163          !                                         ! ================ 
     148         IF(ln_ctl)   CALL prt_ctl(tab2d_1=ui_ice , clinfo1=' lim_dyn  : ui_ice :', tab2d_2=vi_ice , clinfo2=' vi_ice :') 
     149          
     150         ! computation of friction velocity 
     151         ! -------------------------------- 
     152         ! ice-ocean velocity at U & V-points (ui_ice vi_ice at I-point ; ssu_m, ssv_m at U- & V-points) 
     153          
     154         DO jj = 1, jpjm1 
     155            DO ji = 1, fs_jpim1   ! vector opt. 
     156               zu_io(ji,jj) = 0.5 * ( ui_ice(ji+1,jj+1) + ui_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
     157               zv_io(ji,jj) = 0.5 * ( vi_ice(ji+1,jj+1) + vi_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
     158            END DO 
     159         END DO 
     160         ! frictional velocity at T-point 
    164161         DO jj = 2, jpjm1 
    165             zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 
    166             DO ji = 2, jpim1 
    167                ! computation of wind stress over ocean in X and Y direction 
    168 #if defined key_coupled && defined key_lim_cp1 
    169                ztairx =  frld(ji-1,jj  ) * gtaux(ji-1,jj  ) + frld(ji,jj  ) * gtaux(ji,jj  )      & 
    170                   &    + frld(ji-1,jj-1) * gtaux(ji-1,jj-1) + frld(ji,jj-1) * gtaux(ji,jj-1) 
    171  
    172                ztairy =  frld(ji-1,jj  ) * gtauy(ji-1,jj  ) + frld(ji,jj  ) * gtauy(ji,jj  )      & 
    173                   &    + frld(ji-1,jj-1) * gtauy(ji-1,jj-1) + frld(ji,jj-1) * gtauy(ji,jj-1) 
    174 #else 
    175                zsfrld  = frld(ji,jj) + frld(ji-1,jj) + frld(ji-1,jj-1) + frld(ji,jj-1) 
    176                ztairx  = zsfrld * gtaux(ji,jj) 
    177                ztairy  = zsfrld * gtauy(ji,jj) 
    178 #endif 
    179                zsfrldm4 = 4 - frld(ji,jj) - frld(ji-1,jj) - frld(ji-1,jj-1) - frld(ji,jj-1) 
    180                zu_ice   = u_ice(ji,jj) - u_oce(ji,jj) 
    181                zv_ice   = v_ice(ji,jj) - v_oce(ji,jj) 
    182                zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice )  
    183                ztglx   = zsfrldm4 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice )  
    184                ztgly   = zsfrldm4 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice )  
    185  
    186                tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 4 * rau0 ) 
    187                tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 4 * rau0 ) 
     162            DO ji = fs_2, fs_jpim1   ! vector opt. 
     163               ust2s(ji,jj) = 0.5 * cw                                                          & 
     164                  &         * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
     165                  &            + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj) 
    188166            END DO 
    189167         END DO 
    190           
    191          ! computation of friction velocity 
     168         ! 
     169      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
     170         ! 
     171         zcoef = SQRT( 0.5 ) / rau0 
    192172         DO jj = 2, jpjm1 
    193             DO ji = 2, jpim1 
    194  
    195                zu_ice   = u_ice(ji-1,jj-1) - u_oce(ji-1,jj-1) 
    196                zv_ice   = v_ice(ji-1,jj-1) - v_oce(ji-1,jj-1) 
    197                zt11  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
    198  
    199                zu_ice   = u_ice(ji-1,jj) - u_oce(ji-1,jj) 
    200                zv_ice   = v_ice(ji-1,jj) - v_oce(ji-1,jj) 
    201                zt12  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice )  
    202  
    203                zu_ice   = u_ice(ji,jj-1) - u_oce(ji,jj-1) 
    204                zv_ice   = v_ice(ji,jj-1) - v_oce(ji,jj-1) 
    205                zt21  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice )  
    206  
    207                zu_ice   = u_ice(ji,jj) - u_oce(ji,jj) 
    208                zv_ice   = v_ice(ji,jj) - v_oce(ji,jj) 
    209                zt22  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice )  
    210  
    211                ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj) 
    212  
    213                zustm =  ( 1 - frld(ji,jj) ) * 0.25 * ( zt11 + zt12 + zt21 + zt22 )        & 
    214                   &  +        frld(ji,jj)   * SQRT( ztair2 ) 
    215  
    216                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj) 
     173            DO ji = fs_2, fs_jpim1   ! vector opt. 
     174               ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
     175                  &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) 
    217176            END DO 
    218177         END DO 
    219  
    220        ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
    221                      
    222           DO jj = 2, jpjm1 
    223              DO ji = 2, jpim1 
    224 #if defined key_coupled && defined key_lim_cp1 
    225                 tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       & 
    226                    &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 ) 
    227  
    228                 tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       & 
    229                    &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 ) 
    230 #else 
    231                 tio_u(ji,jj) = - gtaux(ji,jj) / rau0 
    232                 tio_v(ji,jj) = - gtauy(ji,jj) / rau0  
    233 #endif 
    234                 ztair2       = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj) 
    235                 zustm        = SQRT( ztair2  ) 
    236  
    237                 ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj) 
    238             END DO 
    239          END DO 
    240  
     178         ! 
    241179      ENDIF 
    242  
     180      ! 
    243181      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point 
    244       CALL lbc_lnk( tio_u, 'I', -1. )   ! I-point (i.e. ice U-V point) 
    245       CALL lbc_lnk( tio_v, 'I', -1. )   ! I-point (i.e. ice U-V point) 
    246  
    247       IF(ln_ctl) THEN  
    248             CALL prt_ctl(tab2d_1=tio_u , clinfo1=' lim_dyn  : tio_u :', tab2d_2=tio_v , clinfo2=' tio_v :') 
    249             CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    250       ENDIF 
     182      ! 
     183      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    251184 
    252185   END SUBROUTINE lim_dyn_2 
     
    257190      !!                  ***  ROUTINE lim_dyn_init_2  *** 
    258191      !! 
    259       !! ** Purpose : Physical constants and parameters linked to the ice 
    260       !!      dynamics 
    261       !! 
    262       !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic 
    263       !!       parameter values called at the first timestep (nit000) 
     192      !! ** Purpose :   Physical constants and parameters linked to the ice 
     193      !!              dynamics 
     194      !! 
     195      !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic 
     196      !!              parameter values 
    264197      !! 
    265198      !! ** input   :   Namelist namicedyn 
    266       !! 
    267       !! history : 
    268       !!  8.5  ! 03-08 (C. Ethe) original code 
    269199      !!------------------------------------------------------------------- 
    270200      NAMELIST/namicedyn/ epsd, alpha,     & 
     
    273203      !!------------------------------------------------------------------- 
    274204 
    275       ! Define the initial parameters 
    276       ! ------------------------- 
    277  
    278       ! Read Namelist namicedyn 
    279       REWIND ( numnam_ice ) 
     205      REWIND ( numnam_ice )                       ! Read Namelist namicedyn 
    280206      READ   ( numnam_ice  , namicedyn ) 
    281       IF(lwp) THEN 
     207 
     208      IF(lwp) THEN                                ! Control print 
    282209         WRITE(numout,*) 
    283210         WRITE(numout,*) 'lim_dyn_init_2: ice parameters for ice dynamics ' 
     
    291218         WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl 
    292219         WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw 
    293          WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg 
     220         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg, ' degrees' 
    294221         WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar 
    295222         WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg 
     
    303230      usecc2 = 1.0 / ( ecc * ecc ) 
    304231      rhoco  = rau0 * cw 
    305       angvg  = angvg * rad 
     232      angvg  = angvg * rad      ! convert angvg from degree to radian 
    306233      sangvg = SIN( angvg ) 
    307234      cangvg = COS( angvg ) 
    308235      pstarh = pstar / 2.0 
    309       sdvt(:,:) = 0.e0 
    310  
    311       !  Diffusion coefficients. 
    312       ahiu(:,:) = ahi0 * umask(:,:,1) 
     236      ! 
     237      ahiu(:,:) = ahi0 * umask(:,:,1)            ! Ice eddy Diffusivity coefficients. 
    313238      ahiv(:,:) = ahi0 * vmask(:,:,1) 
    314  
     239      ! 
    315240   END SUBROUTINE lim_dyn_init_2 
    316241 
  • trunk/NEMO/LIM_SRC_2/limhdf_2.F90

    r823 r888  
    66#if defined key_lim2 
    77   !!---------------------------------------------------------------------- 
    8    !!   'key_lim2'  i                                 LIM 2.0 sea-ice model 
     8   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    99   !!---------------------------------------------------------------------- 
    1010   !!   lim_hdf_2  : diffusion trend on sea-ice variable 
     
    3434   !!---------------------------------------------------------------------- 
    3535   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    36    !! $Header$  
     36   !! $ Id: $ 
    3737   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3838   !!---------------------------------------------------------------------- 
  • trunk/NEMO/LIM_SRC_2/limistate_2.F90

    r823 r888  
    44   !!              Initialisation of diagnostics ice variables 
    55   !!====================================================================== 
    6    !! History :   2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
     6   !! History :   1.0  !  01-04  (C. Ethe, G. Madec)  Original code 
     7   !!             2.0  !  03-08  (G. Madec)  add lim_istate_init 
    78   !!                  !  04-04  (S. Theetten) initialization from a file 
    89   !!                  !  06-07  (S. Masson)  IOM to read the restart 
     10   !!                  !  07-10  (G. Madec)  surface module 
    911   !!-------------------------------------------------------------------- 
    1012#if defined key_lim2 
     
    1820   USE phycst 
    1921   USE ocfzpt 
    20    USE oce             ! dynamics and tracers variables      !!gm used??? 
    21    USE dom_oce                                                     !!gm used??? 
    2222   USE par_ice_2       ! ice parameters 
    2323   USE ice_oce         ! ice variables 
    2424   USE dom_ice_2 
    2525   USE lbclnk 
     26   USE oce 
    2627   USE ice_2 
    2728   USE iom 
     
    4748   !!---------------------------------------------------------------------- 
    4849   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
    49    !! $Header$  
     50   !! $ Id: $ 
    5051   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5152   !!---------------------------------------------------------------------- 
     
    6667      REAL(wp), DIMENSION(jpi,jpj) ::   ztn   ! workspace 
    6768      !-------------------------------------------------------------------- 
    68  
    69        CALL lim_istate_init_2   !  reading the initials parameters of the ice 
    70  
    71       !-- Initialisation of sst,sss,u,v do i=1,jpi 
    72       u_io(:,:)  = 0.e0       ! ice velocity in x direction 
    73       v_io(:,:)  = 0.e0       ! ice velocity in y direction 
    74  
    75       IF( ln_limini ) THEN    !  
    76          
    77          ! Initialisation at tn if no ice or sst_ini if ice 
    78          ! Idem for salinity 
    79  
    80       !--- Criterion for presence (zidto=1.) or absence (zidto=0.) of ice 
    81          DO jj = 1 , jpj 
    82             DO ji = 1 , jpi 
    83                 
    84                zidto = MAX(zzero, - SIGN(1.,frld(ji,jj) - 1.)) 
    85                 
    86                sst_io(ji,jj) = ( nfice - 1 ) * (zidto * sst_ini(ji,jj)  + &   ! use the ocean initial values 
    87                     &          (1.0 - zidto ) * ( tn(ji,jj,1) + rt0 ))        ! tricky trick *(nfice-1) ! 
    88                sss_io(ji,jj) = ( nfice - 1 ) * (zidto * sss_ini(ji,jj) + & 
    89                     &          (1.0 - zidto ) *  sn(ji,jj,1) ) 
    90  
    91                ! to avoid the the melting of ice, several layers (mixed layer) should be 
    92                ! set to sst_ini (sss_ini) if there is ice 
    93                ! example for one layer  
    94                ! tn(ji,jj,1) = zidto * ( sst_ini(ji,jj) - rt0 )  + (1.0 - zidto ) *  tn(ji,jj,1) 
    95                ! sn(ji,jj,1) = zidto * sss_ini(ji,jj)  + (1.0 - zidto ) *  sn(ji,jj,1) 
    96                ! tb(ji,jj,1) = tn(ji,jj,1) 
    97                ! sb(ji,jj,1) = sn(ji,jj,1) 
    98             END DO 
    99          END DO 
    100           
    101           
    102          !  tfu: Melting point of sea water 
    103          tfu(:,:)  = ztf    
    104           
    105          tfu(:,:)  = ABS ( rt0 - 0.0575       * sss_ini(:,:)                               & 
    106               &                    + 1.710523e-03 * sss_ini(:,:) * SQRT( sss_ini(:,:) )    & 
    107               &                    - 2.154996e-04 * sss_ini(:,:) * sss_ini(:,:) ) 
    108       ELSE                     ! 
    109  
     69  
     70      CALL lim_istate_init_2     !  reading the initials parameters of the ice 
     71 
     72      IF( .NOT. ln_limini ) THEN   
    11073          
    11174         ! Initialisation at tn or -2 if ice 
     
    11679            END DO 
    11780         END DO 
    118           
    119          u_io  (:,:) = 0.e0 
    120          v_io  (:,:) = 0.e0 
    121          sst_io(:,:) = ( nfice - 1 ) * ( tn(:,:,1) + rt0 )   ! use the ocean initial values 
    122          sss_io(:,:) = ( nfice - 1 ) *   sn(:,:,1)           ! tricky trick *(nfice-1) ! 
    123           
    124          ! reference salinity 34psu 
     81                   
     82         !  tfu: Melting point of sea water [Kelvin] 
    12583         zs0 = 34.e0 
    126          ztf = ABS ( rt0 - 0.0575       * zs0                           & 
    127               &                    + 1.710523e-03 * zs0 * SQRT( zs0 )   & 
    128               &                    - 2.154996e-04 * zs0 *zs0          ) 
    129           
    130          !  tfu: Melting point of sea water 
    131          tfu(:,:)  = ztf    
     84         ztf = rt0 + ( - 0.0575 + 1.710523e-3 * SQRT( zs0 ) - 2.154996e-4 * zs0 ) * zs0 
     85         tfu(:,:) = ztf 
    13286          
    13387         DO jj = 1, jpj 
     
    152106         tbif  (:,:,2) = tfu(:,:) 
    153107         tbif  (:,:,3) = tfu(:,:) 
    154        
     108 
    155109      ENDIF 
     110      
    156111      fsbbq (:,:)   = 0.e0 
    157112      qstoif(:,:)   = 0.e0 
    158       u_ice (:,:)   = 0.e0 
    159       v_ice (:,:)   = 0.e0 
     113      ui_ice(:,:)   = 0.e0 
     114      vi_ice(:,:)   = 0.e0 
    160115# if defined key_coupled 
    161116      albege(:,:)   = 0.8 * tms(:,:) 
     
    191146 
    192147      CALL lbc_lnk( hsnif, 'T', 1. ) 
    193       CALL lbc_lnk( sist , 'T', 1. ) 
     148      CALL lbc_lnk( sist , 'T', 1. , pval = rt0 )      ! set rt0 on closed boundary (required by bulk formulation) 
    194149      DO jk = 1, jplayersp1 
    195150         CALL lbc_lnk(tbif(:,:,jk), 'T', 1. ) 
     
    197152      CALL lbc_lnk( fsbbq  , 'T', 1. ) 
    198153      CALL lbc_lnk( qstoif , 'T', 1. ) 
    199       CALL lbc_lnk( sss_io , 'T', 1. ) 
    200       ! 
     154 
    201155   END SUBROUTINE lim_istate_2 
    202156 
     
    209163      !! 
    210164      !! ** Method  :   Read the namiceini namelist and check the parameter  
    211       !!                values called at the first timestep (nit000) 
    212       !!                or 
    213       !!                Read 7 variables from a previous restart file 
    214       !!                sst, sst, hicif, hsnif, frld, ts & tbif 
     165      !!       values called at the first timestep (nit000) 
    215166      !! 
    216167      !! ** input   :   Namelist namiceini 
     
    222173         &                hnins, hgins, alins 
    223174      !!------------------------------------------------------------------- 
    224        
    225       ! Read Namelist namiceini  
    226       REWIND ( numnam_ice ) 
     175      ! 
     176      REWIND ( numnam_ice )               ! Read Namelist namiceini  
    227177      READ   ( numnam_ice , namiceini ) 
    228        
    229       IF(.NOT. ln_limini) THEN  
    230          IF(lwp) THEN 
    231             WRITE(numout,*) 
    232             WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 
    233             WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    234             WRITE(numout,*) '         threshold water temp. for initial sea-ice    ttest      = ', ttest 
    235             WRITE(numout,*) '         initial snow thickness in the north          hninn      = ', hninn 
    236             WRITE(numout,*) '         initial ice thickness in the north           hginn      = ', hginn  
    237             WRITE(numout,*) '         initial leads area in the north              alinn      = ', alinn             
    238             WRITE(numout,*) '         initial snow thickness in the south          hnins      = ', hnins  
    239             WRITE(numout,*) '         initial ice thickness in the south           hgins      = ', hgins 
    240             WRITE(numout,*) '         initial leads area in the south              alins      = ', alins 
    241          ENDIF 
     178      ! 
     179      IF(lwp) THEN 
     180         WRITE(numout,*) 
     181         WRITE(numout,*) 'lim_istate_init_2 : ice parameters inititialisation ' 
     182         WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 
     183         WRITE(numout,*) '         threshold water temp. for initial sea-ice    ttest      = ', ttest 
     184         WRITE(numout,*) '         initial snow thickness in the north          hninn      = ', hninn 
     185         WRITE(numout,*) '         initial ice thickness in the north           hginn      = ', hginn  
     186         WRITE(numout,*) '         initial leads area in the north              alinn      = ', alinn             
     187         WRITE(numout,*) '         initial snow thickness in the south          hnins      = ', hnins  
     188         WRITE(numout,*) '         initial ice thickness in the south           hgins      = ', hgins 
     189         WRITE(numout,*) '         initial leads area in the south              alins      = ', alins 
     190         WRITE(numout,*) '         Ice state initialization using input file    ln_limini  = ', ln_limini 
    242191      ENDIF 
    243192 
    244193      IF( ln_limini ) THEN                      ! Ice initialization using input file 
    245  
     194         ! 
    246195         CALL iom_open( 'Ice_initialization.nc', inum_ice ) 
    247  
     196         ! 
    248197         IF( inum_ice > 0 ) THEN 
    249             IF(lwp) THEN 
    250                WRITE(numout,*) ' ' 
    251                WRITE(numout,*) 'lim_istate_init : ice state initialization with : Ice_initialization.nc' 
    252                WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    253                WRITE(numout,*) '         Ice state initialization using input file    ln_limini  = ', ln_limini 
    254                WRITE(numout,*) ' ' 
    255             ENDIF 
     198            IF(lwp) WRITE(numout,*) 
     199            IF(lwp) WRITE(numout,*) '                  ice state initialization with : Ice_initialization.nc' 
    256200             
    257             CALL iom_get( inum_ice, jpdom_data, 'sst'  , sst_ini(:,:) )         
    258             CALL iom_get( inum_ice, jpdom_data, 'sss'  , sss_ini(:,:) )        
    259             CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif  (:,:) )       
    260             CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif  (:,:) )       
    261             CALL iom_get( inum_ice, jpdom_data, 'frld' , frld   (:,:) )      
    262             CALL iom_get( inum_ice, jpdom_data, 'ts'   , sist   (:,:) ) 
     201            CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif )       
     202            CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif )       
     203            CALL iom_get( inum_ice, jpdom_data, 'frld' , frld  )      
     204            CALL iom_get( inum_ice, jpdom_data, 'ts'   , sist  ) 
    263205            CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif(1:nlci,1:nlcj,:),   & 
    264206                 &        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) ) 
     
    268210 
    269211            CALL iom_close( inum_ice) 
    270              
     212            ! 
    271213         ENDIF 
    272214      ENDIF 
    273       ! 
     215      !      
    274216   END SUBROUTINE lim_istate_init_2 
    275217 
  • trunk/NEMO/LIM_SRC_2/limmsh_2.F90

    r823 r888  
    2525   !!---------------------------------------------------------------------- 
    2626   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    27    !! $Header$  
     27   !! $ Id: $ 
    2828   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    2929   !!---------------------------------------------------------------------- 
  • trunk/NEMO/LIM_SRC_2/limrhg_2.F90

    r823 r888  
    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   !!---------------------------------------------------------------------- 
    612#if defined key_lim2 
    713   !!---------------------------------------------------------------------- 
    814   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    915   !!---------------------------------------------------------------------- 
     16   !!---------------------------------------------------------------------- 
    1017   !!   lim_rhg_2   : computes ice velocities 
    1118   !!---------------------------------------------------------------------- 
    12    !! * Modules used 
    13    USE phycst 
    14    USE par_oce 
    15    USE ice_oce         ! ice variables 
    16    USE dom_ice_2 
    17    USE ice_2 
    18    USE lbclnk 
    19    USE lib_mpp 
    20    USE in_out_manager  ! I/O manager 
    21    USE prtctl          ! Print control 
     19   USE par_oce        ! ocean parameter 
     20   USE ice_oce        ! ice variables 
     21   USE sbc_ice        ! surface boundary condition: ice variables 
     22   USE dom_ice_2      ! domaine: ice variables 
     23   USE phycst         ! physical constant 
     24   USE ice_2          ! ice variables 
     25   USE lbclnk         ! lateral boundary condition 
     26   USE lib_mpp        ! MPP library 
     27   USE in_out_manager ! I/O manager 
     28   USE prtctl         ! Print control 
    2229 
    2330   IMPLICIT NONE 
    2431   PRIVATE 
    2532 
    26    !! * Routine accessibility 
    27    PUBLIC lim_rhg_2  ! routine called by lim_dyn_2 
    28  
    29    !! * Module variables 
    30    REAL(wp)  ::           &  ! constant values 
    31       rzero   = 0.e0   ,  & 
    32       rone    = 1.e0 
    33    !!---------------------------------------------------------------------- 
    34    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    35    !! $Header$  
    36    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     33   PUBLIC   lim_rhg_2 ! routine called by lim_dyn 
     34 
     35   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero 
     36   REAL(wp) ::   rone    = 1.e0   !            and  one 
     37 
     38   !! * Substitutions 
     39#  include "vectopt_loop_substitute.h90" 
     40   !!---------------------------------------------------------------------- 
     41   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     42   !! $ Id: $ 
     43   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3744   !!---------------------------------------------------------------------- 
    3845 
     
    4855      !!  viscous-plastic law including shear strength and a bulk rheology. 
    4956      !! 
    50       !! ** Action  : - compute u_ice, v_ice the sea-ice velocity 
     57      !! ** Action  : - compute ui_ice, vi_ice the sea-ice velocity defined 
     58      !!              at I-point 
     59      !!------------------------------------------------------------------- 
     60      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation 
     61      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation 
    5162      !! 
    52       !! History : 
    53       !!   0.0  !  93-12  (M.A. Morales Maqueda.)  Original code 
    54       !!   1.0  !  94-12  (H. Goosse)  
    55       !!   2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp 
     63      INTEGER ::   ji, jj              ! dummy loop indices 
     64      INTEGER ::   iter, jter          ! temporary integers 
     65      CHARACTER (len=50) ::   charout 
     66      REAL(wp) ::   ze11  , ze12  , ze22  , ze21       ! temporary scalars 
     67      REAL(wp) ::   zt11  , zt12  , zt21  , zt22       !    "         " 
     68      REAL(wp) ::   zvis11, zvis21, zvis12, zvis22     !    "         " 
     69      REAL(wp) ::   zgphsx, ztagnx, zunw, zur, zusw    !    "         " 
     70      REAL(wp) ::   zgphsy, ztagny, zvnw, zvr          !    "         " 
     71      REAL(wp) ::   zresm,  za, zac, zmod 
     72      REAL(wp) ::   zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal 
     73      REAL(wp) ::   ztrace2, zdeter, zdelta, zmask, zdgp, zdgi, zdiag 
     74      REAL(wp) ::   za1, zb1, zc1, zd1 
     75      REAL(wp) ::   za2, zb2, zc2, zd2, zden 
     76      REAL(wp) ::   zs11_11, zs11_12, zs11_21, zs11_22 
     77      REAL(wp) ::   zs12_11, zs12_12, zs12_21, zs12_22 
     78      REAL(wp) ::   zs21_11, zs21_12, zs21_21, zs21_22 
     79      REAL(wp) ::   zs22_11, zs22_12, zs22_21, zs22_22 
     80      REAL(wp), DIMENSION(jpi,  jpj  ) ::   zfrld, zmass, zcorl 
     81      REAL(wp), DIMENSION(jpi,  jpj  ) ::   za1ct, za2ct, zresr 
     82      REAL(wp), DIMENSION(jpi,  jpj  ) ::   zc1u, zc1v, zc2u, zc2v 
     83      REAL(wp), DIMENSION(jpi,  jpj  ) ::   zsang 
     84      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zu0, zv0 
     85      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zu_n, zv_n 
     86      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zu_a, zv_a 
     87      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zviszeta, zviseta 
     88      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zzfrld, zztms 
     89      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zi1, zi2, zmasst, zpresh 
     90 
    5691      !!------------------------------------------------------------------- 
    57       ! * Arguments 
    58       INTEGER, INTENT(in) ::   k_j1 ,  &  ! southern j-index for ice computation 
    59          &                     k_jpj      ! northern j-index for ice computation 
    60  
    61       ! * Local variables 
    62       INTEGER ::   ji, jj              ! dummy loop indices 
    63  
    64       INTEGER  :: & 
    65          iim1, ijm1, iip1 , ijp1   , & ! temporary integers 
    66          iter, jter                    !    "          " 
    67  
    68       CHARACTER (len=50) :: charout 
    69  
    70       REAL(wp) :: & 
    71          ze11  , ze12  , ze22  , ze21  ,   &  ! temporary scalars 
    72          zt11  , zt12  , zt21  , zt22  ,   &  !    "         " 
    73          zvis11, zvis21, zvis12, zvis22,   &  !    "         " 
    74          zgphsx, ztagnx, zusw  ,           &  !    "         " 
    75          zgphsy, ztagny                       !    "         " 
    76       REAL(wp) :: & 
    77          zresm, zunw, zvnw, zur, zvr, zmod, za, zac, & 
    78          zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal,  & 
    79          ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag 
    80       REAL(wp),DIMENSION(jpi,jpj) ::   & 
    81          zpresh, zfrld, zmass, zcorl,     & 
    82          zu0, zv0, zviszeta, zviseta,     & 
    83          zc1u, zc1v, zc2u, zc2v, za1ct, za2ct, za1, za2, zb1, zb2,  & 
    84          zc1, zc2, zd1, zd2, zden, zu_ice, zv_ice, zresr 
    85       REAL(wp),DIMENSION(jpi,jpj,2,2) :: & 
    86          zs11, zs12, zs22, zs21 
    87       !!------------------------------------------------------------------- 
     92 
     93!!bug 
     94!!    ui_oce(:,:) = 0.e0 
     95!!    vi_oce(:,:) = 0.e0 
     96!!    write(*,*) 'rhg min, max u & v', maxval(ui_oce), minval(ui_oce), maxval(vi_oce), minval(vi_oce) 
     97!!bug 
    8898       
    8999      !  Store initial velocities 
    90       !  ------------------------ 
    91       zu0(:,:) = u_ice(:,:) 
    92       zv0(:,:) = v_ice(:,:) 
     100      !  ---------------- 
     101      zztms(:,0    ) = 0.e0       ;    zzfrld(:,0    ) = 0.e0 
     102      zztms(:,jpj+1) = 0.e0       ;    zzfrld(:,jpj+1) = 0.e0 
     103      zu0(:,0    ) = 0.e0         ;    zv0(:,0    ) = 0.e0 
     104      zu0(:,jpj+1) = 0.e0         ;    zv0(:,jpj+1) = 0.e0 
     105      zztms(:,1:jpj) = tms(:,:)   ;    zzfrld(:,1:jpj) = frld(:,:) 
     106      zu0(:,1:jpj) = ui_ice(:,:)   ;    zv0(:,1:jpj) = vi_ice(:,:) 
     107 
     108      zu_a(:,:)    = zu0(:,:)     ;   zv_a(:,:) = zv0(:,:) 
     109      zu_n(:,:)    = zu0(:,:)     ;   zv_n(:,:) = zv0(:,:) 
     110 
     111!i 
     112      zi1   (:,:) = 0.e0 
     113      zi2   (:,:) = 0.e0 
     114      zpresh(:,:) = 0.e0 
     115      zmasst(:,:) = 0.e0 
     116!i 
     117!!gm violant 
     118      zfrld(:,:) =0.e0 
     119      zcorl(:,:) =0.e0 
     120      zmass(:,:) =0.e0 
     121      za1ct(:,:) =0.e0 
     122      za2ct(:,:) =0.e0 
     123!!gm end 
     124 
     125      zviszeta(:,:) = 0.e0 
     126      zviseta (:,:) = 0.e0 
     127 
     128!i    zviszeta(:,0    ) = 0.e0    ;    zviseta(:,0    ) = 0.e0 
     129!i    zviszeta(:,jpj  ) = 0.e0    ;    zviseta(:,jpj  ) = 0.e0 
     130!i    zviszeta(:,jpj+1) = 0.e0    ;    zviseta(:,jpj+1) = 0.e0 
     131 
    93132 
    94133      ! Ice mass, ice strength, and wind stress at the center            | 
     
    96135      !------------------------------------------------------------------- 
    97136 
     137!CDIR NOVERRCHK 
    98138      DO jj = k_j1 , k_jpj-1 
     139!CDIR NOVERRCHK 
    99140         DO ji = 1 , jpi 
    100             za1(ji,jj)    = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) 
     141            ! only the sinus changes its sign with the hemisphere 
     142            zsang(ji,jj)  = SIGN( 1.e0, fcor(ji,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere 
     143            ! 
     144            zmasst(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) 
    101145            zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) ) 
    102 #if defined key_lim_cp1 && defined key_coupled 
    103             zb1(ji,jj)    = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    104             zb2(ji,jj)    = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    105 #else 
    106             zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    107             zb2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    108 #endif 
     146!!gm  :: stress given at I-point (F-point for the ocean) only compute the ponderation with the ice fraction (1-frld) 
     147            zi1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 
     148            zi2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 
    109149         END DO 
    110150      END DO 
     
    117157          
    118158      DO jj = k_j1+1, k_jpj-1 
    119          DO ji = 2, jpi 
    120             zstms = tms(ji,jj  ) * wght(ji,jj,2,2) + tms(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    121                &  + tms(ji,jj-1) * wght(ji,jj,2,1) + tms(ji-1,jj-1) * wght(ji,jj,1,1) 
     159         DO ji = fs_2, jpi 
     160            zstms = zztms(ji,jj  ) * wght(ji,jj,2,2) + zztms(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     161               &  + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1) 
    122162            zusw  = 1.0 / MAX( zstms, epsd ) 
    123163 
    124             zt11 = tms(ji  ,jj  ) * frld(ji  ,jj  )  
    125             zt12 = tms(ji-1,jj  ) * frld(ji-1,jj  )  
    126             zt21 = tms(ji  ,jj-1) * frld(ji  ,jj-1)  
    127             zt22 = tms(ji-1,jj-1) * frld(ji-1,jj-1) 
     164            zt11 = zztms(ji  ,jj  ) * zzfrld(ji  ,jj  )  
     165            zt12 = zztms(ji-1,jj  ) * zzfrld(ji-1,jj  )  
     166            zt21 = zztms(ji  ,jj-1) * zzfrld(ji  ,jj-1)  
     167            zt22 = zztms(ji-1,jj-1) * zzfrld(ji-1,jj-1) 
    128168 
    129169            ! Leads area. 
     
    131171               &             + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw 
    132172 
    133             ! Mass and coriolis coeff. 
    134             zmass(ji,jj) = ( za1(ji,jj  ) * wght(ji,jj,2,2) + za1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    135                &           + za1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
     173            ! Mass and coriolis coeff. at I-point 
     174            zmass(ji,jj) = ( zmasst(ji,jj  ) * wght(ji,jj,2,2) + zmasst(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     175               &           + zmasst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
    136176            zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj) 
    137177 
    138178            ! Wind stress. 
    139 #if defined key_lim_cp1 && defined key_coupled 
    140             ztagnx = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    141                &     + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
    142             ztagny = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    143                &     + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
    144 #else 
    145             ztagnx = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    146                &     + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj) 
    147             ztagny = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    148                &     + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj) 
    149 #endif 
     179            ! always provide stress at I-point (ocean F-point) 
     180            ztagnx = ( zi1(ji,jj  ) * wght(ji,jj,2,2) + zi1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     181               &     + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utaui_ice(ji,jj) 
     182            ztagny = ( zi2(ji,jj  ) * wght(ji,jj,2,2) + zi2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     183               &     + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtaui_ice(ji,jj) 
    150184 
    151185            ! Gradient of ice strength 
     
    161195 
    162196            ! Computation of the velocity field taking into account the ice-ice interaction.                                  
    163             ! Terms that are independent of the velocity field. 
    164             za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_oce(ji,jj) - zgphsx 
    165             za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_oce(ji,jj) - zgphsy 
     197            ! Terms that are independent of the ice velocity field. 
     198            za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * vi_oce(ji,jj) - zgphsx 
     199            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * ui_oce(ji,jj) - zgphsy 
    166200         END DO 
    167201      END DO 
    168  
    169 !! inutile!! 
    170 !!??    CALL lbc_lnk( za1ct, 'I', -1. ) 
    171 !!??    CALL lbc_lnk( za2ct, 'I', -1. ) 
    172202 
    173203 
     
    182212         ! Computation of free drift field for free slip boundary conditions. 
    183213 
    184            DO jj = k_j1, k_jpj-1 
    185               DO ji = 1, jpim1 
    186                  !- Rate of strain tensor. 
    187                  zt11 =   akappa(ji,jj,1,1) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) - u_ice(ji,jj  ) - u_ice(ji  ,jj+1) )  & 
    188                     &   + akappa(ji,jj,1,2) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) + v_ice(ji,jj  ) + v_ice(ji  ,jj+1) ) 
    189                  zt12 = - akappa(ji,jj,2,2) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) - u_ice(ji,jj+1) - u_ice(ji+1,jj+1) )  & 
    190                     &   - akappa(ji,jj,2,1) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) + v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) 
    191                  zt22 = - akappa(ji,jj,2,2) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) - v_ice(ji,jj+1) - v_ice(ji+1,jj+1) )  & 
    192                     &   + akappa(ji,jj,2,1) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) + u_ice(ji,jj+1) + u_ice(ji+1,jj+1) ) 
    193                  zt21 =   akappa(ji,jj,1,1) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) - v_ice(ji,jj  ) - v_ice(ji  ,jj+1) )  & 
    194                     &   - akappa(ji,jj,1,2) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) + u_ice(ji,jj  ) + u_ice(ji  ,jj+1) ) 
    195  
    196                  !- Rate of strain tensor.  
    197                  zdgp = zt11 + zt22 
    198                  zdgi = zt12 + zt21 
    199                  ztrace2 = zdgp * zdgp  
    200                  zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi 
    201  
    202                  !  Creep limit depends on the size of the grid. 
    203                  zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2),  creepl) 
    204  
    205                  !-  Computation of viscosities. 
    206                  zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 
    207                  zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 
    208               END DO 
    209            END DO 
    210 !!??       CALL lbc_lnk( zviszeta, 'I', -1. )  ! or T point???   semble reellement inutile 
    211 !!??       CALL lbc_lnk( zviseta , 'I', -1. ) 
    212  
    213  
    214            !-  Determination of zc1u, zc2u, zc1v and zc2v. 
    215            DO jj = k_j1+1, k_jpj-1 
    216               DO ji = 2, jpim1 
    217                  ze11   =  akappa(ji-1,jj-1,1,1) 
    218                  ze12   = +akappa(ji-1,jj-1,2,2) 
    219                  ze22   =  akappa(ji-1,jj-1,2,1) 
    220                  ze21   = -akappa(ji-1,jj-1,1,2) 
    221                  zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
    222                  zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
    223                  zvis12 =       zviseta (ji-1,jj-1) + dm 
    224                  zvis21 =       zviseta (ji-1,jj-1) 
    225  
    226                  zdiag = zvis22 * ( ze11 + ze22 ) 
    227                  zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag 
    228                  zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21 
    229                  zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag 
    230                  zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12 
    231  
    232                  ze11   = -akappa(ji,jj-1,1,1) 
    233                  ze12   = +akappa(ji,jj-1,2,2) 
    234                  ze22   =  akappa(ji,jj-1,2,1) 
    235                  ze21   = -akappa(ji,jj-1,1,2) 
    236                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
    237                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
    238                  zvis12 =       zviseta (ji,jj-1) + dm 
    239                  zvis21 =       zviseta (ji,jj-1) 
    240  
    241                  zdiag = zvis22 * ( ze11 + ze22 ) 
    242                  zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag 
    243                  zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21 
    244                  zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag 
    245                  zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12 
    246  
    247                  ze11   =  akappa(ji-1,jj,1,1) 
    248                  ze12   = -akappa(ji-1,jj,2,2) 
    249                  ze22   =  akappa(ji-1,jj,2,1) 
    250                  ze21   = -akappa(ji-1,jj,1,2) 
    251                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
    252                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
    253                  zvis12 =       zviseta (ji-1,jj) + dm 
    254                  zvis21 =       zviseta (ji-1,jj) 
    255  
    256                  zdiag = zvis22 * ( ze11 + ze22 )  
    257                  zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag 
    258                  zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21 
    259                  zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag 
    260                  zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12 
    261  
    262                  ze11   = -akappa(ji,jj,1,1) 
    263                  ze12   = -akappa(ji,jj,2,2) 
    264                  ze22   =  akappa(ji,jj,2,1) 
    265                  ze21   = -akappa(ji,jj,1,2) 
    266                  zvis11 = 2.0 * zviseta (ji,jj) + dm 
    267                  zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
    268                  zvis12 =       zviseta (ji,jj) + dm 
    269                  zvis21 =       zviseta (ji,jj) 
    270  
    271                  zdiag = zvis22 * ( ze11 + ze22 ) 
    272                  zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag 
    273                  zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21 
    274                  zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag  
    275                  zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12 
    276               END DO 
    277            END DO 
    278  
    279            DO jj = k_j1+1, k_jpj-1 
    280               DO ji = 2, jpim1 
    281                  zc1u(ji,jj) =   & 
    282                     + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   & 
    283                     - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   & 
    284                     - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   & 
    285                     + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   & 
    286                     + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   & 
    287                     + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   & 
    288                     - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   & 
    289                     - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 
    290                   
    291                  zc2u(ji,jj) =   & 
    292                     + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   & 
    293                     - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   & 
    294                     - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   & 
    295                     + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   & 
    296                     - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   & 
    297                     - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   & 
    298                     + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   & 
    299                     + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 
    300              END DO 
    301            END DO 
    302  
    303            DO jj = k_j1+1, k_jpj-1 
    304               DO ji = 2, jpim1 
    305                  !  zc1v , zc2v. 
    306                  ze11   =  akappa(ji-1,jj-1,1,2) 
    307                  ze12   = -akappa(ji-1,jj-1,2,1) 
    308                  ze22   = +akappa(ji-1,jj-1,2,2) 
    309                  ze21   =  akappa(ji-1,jj-1,1,1) 
    310                  zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
    311                  zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
    312                  zvis12 =       zviseta (ji-1,jj-1) + dm 
    313                  zvis21 =       zviseta (ji-1,jj-1) 
    314  
    315                  zdiag = zvis22 * ( ze11 + ze22 ) 
    316                  zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag 
    317                  zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21 
    318                  zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag 
    319                  zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12 
     214!CDIR NOVERRCHK 
     215         DO jj = k_j1, k_jpj-1 
     216!CDIR NOVERRCHK 
     217            DO ji = 1, fs_jpim1 
     218               !- Rate of strain tensor. 
     219               zt11 =   akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji,jj  ) - zu_a(ji  ,jj+1) )  & 
     220                  &   + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji,jj  ) + zv_a(ji  ,jj+1) ) 
     221               zt12 = - akappa(ji,jj,2,2) * ( zu_a(ji  ,jj) + zu_a(ji+1,jj  ) - zu_a(ji,jj+1) - zu_a(ji+1,jj+1) )  & 
     222                  &   - akappa(ji,jj,2,1) * ( zv_a(ji  ,jj) + zv_a(ji+1,jj  ) + zv_a(ji,jj+1) + zv_a(ji+1,jj+1) ) 
     223               zt22 = - akappa(ji,jj,2,2) * ( zv_a(ji  ,jj) + zv_a(ji+1,jj  ) - zv_a(ji,jj+1) - zv_a(ji+1,jj+1) )  & 
     224                  &   + akappa(ji,jj,2,1) * ( zu_a(ji  ,jj) + zu_a(ji+1,jj  ) + zu_a(ji,jj+1) + zu_a(ji+1,jj+1) ) 
     225               zt21 =   akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji,jj  ) - zv_a(ji  ,jj+1) )  & 
     226                  &   - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji,jj  ) + zu_a(ji  ,jj+1) ) 
     227 
     228               !- Rate of strain tensor.  
     229               zdgp = zt11 + zt22 
     230               zdgi = zt12 + zt21 
     231               ztrace2 = zdgp * zdgp  
     232               zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi 
     233 
     234               !  Creep limit depends on the size of the grid. 
     235               zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2 ),  creepl) 
     236 
     237               !-  Computation of viscosities. 
     238               zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 
     239               zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 
     240            END DO 
     241         END DO 
     242 
     243         !-  Determination of zc1u, zc2u, zc1v and zc2v. 
     244         DO jj = k_j1+1, k_jpj-1 
     245            DO ji = fs_2, fs_jpim1 
     246               !* zc1u , zc2v 
     247               zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
     248               zvis12 =       zviseta (ji-1,jj-1) + dm 
     249               zvis21 =       zviseta (ji-1,jj-1) 
     250               zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
     251               zdiag  = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) ) 
     252               zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,1) + zdiag 
     253               zs12_11 =  zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2) 
     254               zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2) 
     255               zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,1) + zdiag 
     256 
     257               zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
     258               zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
     259               zvis12 =       zviseta (ji,jj-1) + dm 
     260               zvis21 =       zviseta (ji,jj-1) 
     261               zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) ) 
     262               zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag 
     263               zs12_21 =  zvis12 * akappa(ji,jj-1,2,2) - zvis21 * akappa(ji,jj-1,1,2) 
     264               zs22_21 =  zvis11 * akappa(ji,jj-1,2,1) + zdiag 
     265               zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2) 
     266 
     267               zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
     268               zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
     269               zvis12 =       zviseta (ji-1,jj) + dm 
     270               zvis21 =       zviseta (ji-1,jj) 
     271               zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) ) 
     272               zs11_12 =  zvis11 * akappa(ji-1,jj,1,1) + zdiag 
     273               zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2) 
     274               zs22_12 =  zvis11 * akappa(ji-1,jj,2,1) + zdiag 
     275               zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2) 
     276 
     277               zvis11 = 2.0 * zviseta (ji,jj) + dm 
     278               zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
     279               zvis12 =       zviseta (ji,jj) + dm 
     280               zvis21 =       zviseta (ji,jj) 
     281               zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) ) 
     282               zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag 
     283               zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2) 
     284               zs22_22 =  zvis11 * akappa(ji,jj,2,1) + zdiag 
     285               zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2) 
     286 
     287               zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22   & 
     288                  &          - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12   & 
     289                  &          - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11   & 
     290                  &          + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12   & 
     291                  &          + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21   & 
     292                  &          + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22   & 
     293                  &          - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21   & 
     294                  &          - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 
     295 
     296               zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22   & 
     297                  &          - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12   & 
     298                  &          - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11   & 
     299                  &          + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12   & 
     300                  &          - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21   & 
     301                  &          - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22   & 
     302                  &          + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21   & 
     303                  &          + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 
     304 
     305               !* zc1v , zc2v. 
     306               zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
     307               zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
     308               zvis12 =       zviseta (ji-1,jj-1) + dm 
     309               zvis21 =       zviseta (ji-1,jj-1) 
     310               zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) ) 
     311               zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,2) + zdiag 
     312               zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1) 
     313               zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,2) + zdiag 
     314               zs21_11 =  zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1) 
    320315  
    321                  ze11   =  akappa(ji,jj-1,1,2) 
    322                  ze12   = -akappa(ji,jj-1,2,1) 
    323                  ze22   = +akappa(ji,jj-1,2,2) 
    324                  ze21   = -akappa(ji,jj-1,1,1) 
    325                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
    326                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
    327                  zvis12 =       zviseta (ji,jj-1) + dm 
    328                  zvis21 =       zviseta (ji,jj-1) 
    329  
    330                  zdiag = zvis22 * ( ze11 + ze22 ) 
    331                  zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag 
    332                  zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21 
    333                  zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag 
    334                  zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12 
    335  
    336                  ze11   =  akappa(ji-1,jj,1,2) 
    337                  ze12   = -akappa(ji-1,jj,2,1) 
    338                  ze22   = -akappa(ji-1,jj,2,2) 
    339                  ze21   =  akappa(ji-1,jj,1,1) 
    340                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
    341                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
    342                  zvis12 =       zviseta (ji-1,jj) + dm 
    343                  zvis21 =       zviseta (ji-1,jj) 
    344  
    345                  zdiag = zvis22 * ( ze11 + ze22 ) 
    346                  zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag 
    347                  zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21 
    348                  zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag 
    349                  zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12 
    350  
    351                  ze11   =  akappa(ji,jj,1,2) 
    352                  ze12   = -akappa(ji,jj,2,1) 
    353                  ze22   = -akappa(ji,jj,2,2) 
    354                  ze21   = -akappa(ji,jj,1,1) 
    355                  zvis11 = 2.0 * zviseta (ji,jj) + dm 
    356                  zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
    357                  zvis12 =       zviseta (ji,jj) + dm 
    358                  zvis21 =       zviseta (ji,jj) 
    359  
    360                  zdiag = zvis22 * ( ze11 + ze22 ) 
    361                  zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag 
    362                  zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21 
    363                  zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag 
    364                  zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12 
    365  
    366               END DO 
    367            END DO 
    368  
    369            DO jj = k_j1+1, k_jpj-1 
    370               DO ji = 2, jpim1 
    371                  zc1v(ji,jj) =   & 
    372                     + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   & 
    373                     - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   & 
    374                     - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   & 
    375                     + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   & 
    376                     + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   & 
    377                     + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   & 
    378                     - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   & 
    379                     - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 
    380                  zc2v(ji,jj) =   & 
    381                     + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   & 
    382                     - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   & 
    383                     - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   & 
    384                     + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   & 
    385                     - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   & 
    386                     - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   & 
    387                     + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   & 
    388                     + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 
    389               END DO 
    390            END DO 
    391  
    392          ! Relaxation. 
    393             
    394 iflag:   DO jter = 1 , nbitdr 
    395  
    396             !  Store previous drift field.    
    397             DO jj = k_j1, k_jpj-1 
    398                zu_ice(:,jj) = u_ice(:,jj) 
    399                zv_ice(:,jj) = v_ice(:,jj) 
     316               zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
     317               zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
     318               zvis12 =       zviseta (ji,jj-1) + dm 
     319               zvis21 =       zviseta (ji,jj-1) 
     320               zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) ) 
     321               zs11_21 =  zvis11 * akappa(ji,jj-1,1,2) + zdiag 
     322               zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1) 
     323               zs22_21 =  zvis11 * akappa(ji,jj-1,2,2) + zdiag 
     324               zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1) 
     325 
     326               zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
     327               zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
     328               zvis12 =       zviseta (ji-1,jj) + dm 
     329               zvis21 =       zviseta (ji-1,jj) 
     330               zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) ) 
     331               zs11_12 =  zvis11 * akappa(ji-1,jj,1,2) + zdiag 
     332               zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1) 
     333               zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag 
     334               zs21_12 =  zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1) 
     335 
     336               zvis11 = 2.0 * zviseta (ji,jj) + dm 
     337               zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
     338               zvis12 =       zviseta (ji,jj) + dm 
     339               zvis21 =       zviseta (ji,jj) 
     340               zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) ) 
     341               zs11_22 =  zvis11 * akappa(ji,jj,1,2) + zdiag 
     342               zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1) 
     343               zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag 
     344               zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1) 
     345 
     346               zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22   & 
     347                  &          - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12   & 
     348                  &          - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11   & 
     349                  &          + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12   & 
     350                  &          + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21   & 
     351                  &          + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22   & 
     352                  &          - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21   & 
     353                  &          - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 
     354 
     355               zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22   & 
     356                  &          - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12   & 
     357                  &          - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11   & 
     358                  &          + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12   & 
     359                  &          - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21   & 
     360                  &          - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22   & 
     361                  &          + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21   & 
     362                  &          + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 
    400363            END DO 
    401  
     364         END DO 
     365 
     366         ! GAUSS-SEIDEL method 
     367         !                                                      ! ================ ! 
     368iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! 
     369            !                                                   ! ================ ! 
     370!CDIR NOVERRCHK 
    402371            DO jj = k_j1+1, k_jpj-1 
    403                zsang  = SIGN( 1.e0, fcor(1,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere 
    404                DO ji = 2, jpim1 
    405                  zur     = u_ice(ji,jj) - u_oce(ji,jj) 
    406                  zvr     = v_ice(ji,jj) - v_oce(ji,jj) 
    407                  zmod    = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) ) 
    408                  za      = rhoco * zmod 
    409                  zac     = za * cangvg 
    410                   zmpzas  = alpha * zcorl(ji,jj) + za * zsang 
     372!CDIR NOVERRCHK 
     373               DO ji = fs_2, fs_jpim1 
     374                  ! 
     375                  ze11 =   akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj) 
     376                  ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj) 
     377                  ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj) 
     378                  ze21 =   akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj) 
     379                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
     380                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
     381                  zvis12 =       zviseta (ji,jj-1) + dm 
     382                  zvis21 =       zviseta (ji,jj-1) 
     383                  zdiag = zvis22 * ( ze11 + ze22 ) 
     384                  zs11_21 =  zvis11 * ze11 + zdiag 
     385                  zs12_21 =  zvis12 * ze12 + zvis21 * ze21 
     386                  zs22_21 =  zvis11 * ze22 + zdiag 
     387                  zs21_21 =  zvis12 * ze21 + zvis21 * ze12 
     388 
     389                  ze11 =   akappa(ji-1,jj,1,1) * ( zu_a(ji  ,jj+1) - zu_a(ji-1,jj+1) )   & 
     390                     &   + akappa(ji-1,jj,1,2) * ( zv_a(ji  ,jj+1) + zv_a(ji-1,jj+1) ) 
     391                  ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji  ,jj+1) )   & 
     392                     &   - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji  ,jj+1) ) 
     393                  ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji  ,jj+1) )   & 
     394                     &   + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji  ,jj+1) ) 
     395                  ze21 =   akappa(ji-1,jj,1,1) * ( zv_a(ji  ,jj+1) - zv_a(ji-1,jj+1) )   & 
     396                     &   - akappa(ji-1,jj,1,2) * ( zu_a(ji  ,jj+1) + zu_a(ji-1,jj+1) ) 
     397                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
     398                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
     399                  zvis12 =       zviseta (ji-1,jj) + dm 
     400                  zvis21 =       zviseta (ji-1,jj) 
     401                  zdiag = zvis22 * ( ze11 + ze22 ) 
     402                  zs11_12 =  zvis11 * ze11 + zdiag 
     403                  zs12_12 =  zvis12 * ze12 + zvis21 * ze21 
     404                  zs22_12 =  zvis11 * ze22 + zdiag 
     405                  zs21_12 =  zvis12 * ze21 + zvis21 * ze12 
     406 
     407                  ze11 =   akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji  ,jj+1) )   & 
     408                     &   + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji  ,jj+1) ) 
     409                  ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji  ,jj+1) - zu_a(ji+1,jj+1) )   & 
     410                     &   - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji  ,jj+1) + zv_a(ji+1,jj+1) ) 
     411                  ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji  ,jj+1) - zv_a(ji+1,jj+1) )   & 
     412                     &   + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji  ,jj+1) + zu_a(ji+1,jj+1) ) 
     413                  ze21 =   akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji  ,jj+1) )   & 
     414                     &   - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji  ,jj+1) ) 
     415                  zvis11 = 2.0 * zviseta (ji,jj) + dm 
     416                  zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
     417                  zvis12 =       zviseta (ji,jj) + dm 
     418                  zvis21 =       zviseta (ji,jj) 
     419                  zdiag = zvis22 * ( ze11 + ze22 ) 
     420                  zs11_22 =  zvis11 * ze11 + zdiag 
     421                  zs12_22 =  zvis12 * ze12 + zvis21 * ze21 
     422                  zs22_22 =  zvis11 * ze22 + zdiag 
     423                  zs21_22 =  zvis12 * ze21 + zvis21 * ze12 
     424 
     425            ! 2nd part 
     426                  ze11 =   akappa(ji-1,jj-1,1,1) * ( zu_a(ji  ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) )   & 
     427                     &   + akappa(ji-1,jj-1,1,2) * ( zv_a(ji  ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) ) 
     428                  ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji  ,jj-1) - zu_a(ji-1,jj) )   & 
     429                     &   - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji  ,jj-1) + zv_a(ji-1,jj) ) 
     430                  ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji  ,jj-1) - zv_a(ji-1,jj) )   & 
     431                     &   + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji  ,jj-1) + zu_a(ji-1,jj) ) 
     432                  ze21 =   akappa(ji-1,jj-1,1,1) * ( zv_a(ji  ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) )   & 
     433                     &  -  akappa(ji-1,jj-1,1,2) * ( zu_a(ji  ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) ) 
     434                  zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
     435                  zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
     436                  zvis12 =       zviseta (ji-1,jj-1) + dm 
     437                  zvis21 =       zviseta (ji-1,jj-1) 
     438                  zdiag = zvis22 * ( ze11 + ze22 ) 
     439                  zs11_11 =  zvis11 * ze11 + zdiag 
     440                  zs12_11 =  zvis12 * ze12 + zvis21 * ze21 
     441                  zs22_11 =  zvis11 * ze22 + zdiag 
     442                  zs21_11 =  zvis12 * ze21 + zvis21 * ze12 
     443 
     444                  ze11 =   akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji  ,jj-1) )   & 
     445                     &   + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji  ,jj-1) ) 
     446                  ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji  ,jj-1) + zu_a(ji+1,jj-1) )   & 
     447                     &   - akappa(ji,jj-1,2,1) * ( zv_a(ji  ,jj-1) + zv_a(ji+1,jj-1) ) 
     448                  ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji  ,jj-1) + zv_a(ji+1,jj-1) )   & 
     449                     &   + akappa(ji,jj-1,2,1) * ( zu_a(ji  ,jj-1) + zu_a(ji+1,jj-1) ) 
     450                  ze21 =   akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji  ,jj-1) )   & 
     451                     &   - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji  ,jj-1) ) 
     452                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
     453                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
     454                  zvis12 =       zviseta (ji,jj-1) + dm 
     455                  zvis21 =       zviseta (ji,jj-1) 
     456                  zdiag = zvis22 * ( ze11 + ze22 ) 
     457                  zs11_21 =  zs11_21 + zvis11 * ze11 + zdiag 
     458                  zs12_21 =  zs12_21 + zvis12 * ze12 + zvis21 * ze21 
     459                  zs22_21 =  zs22_21 + zvis11 * ze22 + zdiag 
     460                  zs21_21 =  zs21_21 + zvis12 * ze21 + zvis21 * ze12 
     461 
     462                  ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj) 
     463                  ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj) 
     464                  ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj) 
     465                  ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj) 
     466                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
     467                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
     468                  zvis12 =       zviseta (ji-1,jj) + dm 
     469                  zvis21 =       zviseta (ji-1,jj) 
     470                  zdiag = zvis22 * ( ze11 + ze22 ) 
     471                  zs11_12 =  zs11_12 + zvis11 * ze11 + zdiag 
     472                  zs12_12 =  zs12_12 + zvis12 * ze12 + zvis21 * ze21 
     473                  zs22_12 =  zs22_12 + zvis11 * ze22 + zdiag 
     474                  zs21_12 =  zs21_12 + zvis12 * ze21 + zvis21 * ze12 
     475 
     476                  zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22  & 
     477                     &  - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12  & 
     478                     &  - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11  & 
     479                     &  + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12  & 
     480                     &  + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21  & 
     481                     &  + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22  & 
     482                     &  - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21  & 
     483                     &  - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 
     484 
     485                  zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22  & 
     486                     &  - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12  & 
     487                     &  - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11  & 
     488                     &  + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12  & 
     489                     &  - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21  & 
     490                     &  - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22  & 
     491                     &  + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21  & 
     492                     &  + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 
     493 
     494                  zur     = zu_a(ji,jj) - ui_oce(ji,jj) 
     495                  zvr     = zv_a(ji,jj) - vi_oce(ji,jj) 
     496!!!! 
     497                  zmod    = SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 
     498                  za      = rhoco * zmod 
     499!!!! 
     500!!gm chg resul    za      = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 
     501                  zac     = za * cangvg 
     502                  zmpzas  = alpha * zcorl(ji,jj) + za * zsang(ji,jj) 
    411503                  zmassdt = zusdtp * zmass(ji,jj) 
    412504                  zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj) 
    413505 
    414                   za1(ji,jj) =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   & 
    415                      &        + za * ( cangvg * u_oce(ji,jj) - zsang * v_oce(ji,jj) ) 
    416  
    417                   za2(ji,jj) =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   & 
    418                      &        + za * ( cangvg * v_oce(ji,jj) + zsang * u_oce(ji,jj) ) 
    419  
    420                   zb1(ji,jj)  = zmassdt + zac - zc1u(ji,jj) 
    421                   zb2(ji,jj)  = zmpzas  - zc2u(ji,jj) 
    422                   zc1(ji,jj)  = zmpzas  + zc1v(ji,jj) 
    423                   zc2(ji,jj)  = zmassdt + zac  - zc2v(ji,jj)  
    424                   zdeter      = zc1(ji,jj) * zb2(ji,jj) + zc2(ji,jj) * zb1(ji,jj) 
    425                   zden(ji,jj) = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 
     506                  za1 =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   & 
     507                     &        + za * ( cangvg * ui_oce(ji,jj) - zsang(ji,jj) * vi_oce(ji,jj) ) 
     508                  za2 =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   & 
     509                     &        + za * ( cangvg * vi_oce(ji,jj) + zsang(ji,jj) * ui_oce(ji,jj) ) 
     510                  zb1    = zmassdt + zac - zc1u(ji,jj) 
     511                  zb2    = zmpzas        - zc2u(ji,jj) 
     512                  zc1    = zmpzas        + zc1v(ji,jj) 
     513                  zc2    = zmassdt + zac - zc2v(ji,jj) 
     514                  zdeter = zc1 * zb2 + zc2 * zb1 
     515                  zden   = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 
     516                  zunw   = (  ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden 
     517                  zvnw   = (  ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden 
     518                  zmask  = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 
     519 
     520                  zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask 
     521                  zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask 
    426522               END DO 
    427523            END DO 
    428524 
    429             ! The computation of ice interaction term is splitted into two parts 
    430             !------------------------------------------------------------------------- 
    431  
    432             ! Terms that do not involve already up-dated velocities. 
    433           
    434             DO jj = k_j1+1, k_jpj-1 
    435                DO ji = 2, jpim1 
    436                   iim1 = ji 
    437                   ijm1 = jj - 1 
    438                   iip1 = ji + 1 
    439                   ijp1 = jj 
    440                   ze11 =   akappa(iim1,ijm1,1,1) * u_ice(iip1,ijp1) + akappa(iim1,ijm1,1,2) * v_ice(iip1,ijp1) 
    441                   ze12 = + akappa(iim1,ijm1,2,2) * u_ice(iip1,ijp1) - akappa(iim1,ijm1,2,1) * v_ice(iip1,ijp1) 
    442                   ze22 = + akappa(iim1,ijm1,2,2) * v_ice(iip1,ijp1) + akappa(iim1,ijm1,2,1) * u_ice(iip1,ijp1) 
    443                   ze21 =   akappa(iim1,ijm1,1,1) * v_ice(iip1,ijp1) - akappa(iim1,ijm1,1,2) * u_ice(iip1,ijp1) 
    444                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    445                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    446                   zvis12 =       zviseta (iim1,ijm1) + dm 
    447                   zvis21 =       zviseta (iim1,ijm1) 
    448                   zdiag = zvis22 * ( ze11 + ze22 ) 
    449                   zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag 
    450                   zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21 
    451                   zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag 
    452                   zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12 
    453  
    454  
    455                   iim1 = ji - 1 
    456                   ijm1 = jj 
    457                   iip1 = ji 
    458                   ijp1 = jj + 1                    
    459                   ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   & 
    460                      &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijp1) + v_ice(iim1,ijp1) ) 
    461                   ze12 = + akappa(iim1,ijm1,2,2) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )   & 
    462                      &   - akappa(iim1,ijm1,2,1) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) 
    463                   ze22 = + akappa(iim1,ijm1,2,2) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )   & 
    464                      &   + akappa(iim1,ijm1,2,1) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) 
    465                   ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   & 
    466                      &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 
    467                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    468                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    469                   zvis12 =       zviseta (iim1,ijm1) + dm 
    470                   zvis21 =       zviseta (iim1,ijm1) 
    471                   zdiag = zvis22 * ( ze11 + ze22 ) 
    472                   zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag 
    473                   zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21 
    474                   zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag 
    475                   zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12 
    476  
    477                   iim1 = ji 
    478                   ijm1 = jj 
    479                   iip1 = ji + 1 
    480                   ijp1 = jj + 1 
    481                   ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   & 
    482                      &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) + v_ice(iim1,ijp1) ) 
    483                   ze12 = - akappa(iim1,ijm1,2,2) * ( u_ice(iip1,ijm1) - u_ice(iim1,ijp1) - u_ice(iip1,ijp1) )   & 
    484                      &   - akappa(iim1,ijm1,2,1) * ( v_ice(iip1,ijm1) + v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) 
    485                   ze22 = - akappa(iim1,ijm1,2,2) * ( v_ice(iip1,ijm1) - v_ice(iim1,ijp1) - v_ice(iip1,ijp1) )   & 
    486                      &   + akappa(iim1,ijm1,2,1) * ( u_ice(iip1,ijm1) + u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) 
    487                   ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   & 
    488                      &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) + u_ice(iim1,ijp1) )  
    489                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    490                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    491                   zvis12 =       zviseta (iim1,ijm1) + dm 
    492                   zvis21 =       zviseta (iim1,ijm1) 
    493  
    494                   zdiag = zvis22 * ( ze11 + ze22 ) 
    495                   zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag 
    496                   zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21 
    497                   zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag 
    498                   zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12 
    499  
    500                END DO 
     525            CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. ) 
     526            CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. ) 
     527 
     528            ! Test of Convergence 
     529            DO jj = k_j1+1 , k_jpj-1 
     530               zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) ) 
    501531            END DO 
    502  
    503             ! Terms involving already up-dated velocities. 
    504             !-Using the arrays zu_ice and zv_ice in the computation of the terms ze leads to JACOBI's method;  
    505             ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method. 
    506               
    507             DO jj = k_j1+1, k_jpj-1 
    508                DO ji = 2, jpim1 
    509                   iim1 = ji - 1 
    510                   ijm1 = jj - 1 
    511                   iip1 = ji 
    512                   ijp1 = jj 
    513                   ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) - zu_ice(iim1,ijp1) )   & 
    514                      &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) + zv_ice(iim1,ijp1) ) 
    515                   ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) - zu_ice(iim1,ijp1) )   & 
    516                      &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) + zv_ice(iim1,ijp1) ) 
    517                   ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) - zv_ice(iim1,ijp1) )   & 
    518                      &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) + zu_ice(iim1,ijp1) ) 
    519                   ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) - zv_ice(iim1,ijp1) )   & 
    520                      &  -  akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) + zu_ice(iim1,ijp1) ) 
    521                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    522                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    523                   zvis12 =       zviseta (iim1,ijm1) + dm 
    524                   zvis21 =       zviseta (iim1,ijm1) 
    525  
    526                   zdiag = zvis22 * ( ze11 + ze22 ) 
    527                   zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag 
    528                   zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21 
    529                   zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag 
    530                   zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12 
    531  
    532 #if defined key_agrif 
    533              END DO 
    534           END DO 
    535  
    536           DO jj = k_j1+1, k_jpj-1 
    537              DO ji = 2, jpim1 
    538 #endif 
    539  
    540                   iim1 = ji 
    541                   ijm1 = jj - 1 
    542                   iip1 = ji + 1 
    543                   ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) )   & 
    544                      &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) ) 
    545                   ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )   & 
    546                      &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) ) 
    547                   ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )   & 
    548                      &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) ) 
    549                   ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) )   & 
    550                      &   - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) ) 
    551                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    552                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    553                   zvis12 =       zviseta (iim1,ijm1) + dm 
    554                   zvis21 =       zviseta (iim1,ijm1) 
    555  
    556                   zdiag = zvis22 * ( ze11 + ze22 ) 
    557                   zs11(ji,jj,2,1) =  zs11(ji,jj,2,1) + zvis11 * ze11 + zdiag 
    558                   zs12(ji,jj,2,1) =  zs12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21 
    559                   zs22(ji,jj,2,1) =  zs22(ji,jj,2,1) + zvis11 * ze22 + zdiag 
    560                   zs21(ji,jj,2,1) =  zs21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12 
    561  
    562  
    563                   iim1 = ji - 1 
    564                   ijm1 = jj  
    565                   ze11 = - akappa(iim1,ijm1,1,1) * zu_ice(iim1,ijm1) + akappa(iim1,ijm1,1,2) * zv_ice(iim1,ijm1) 
    566                   ze12 = - akappa(iim1,ijm1,2,2) * zu_ice(iim1,ijm1) - akappa(iim1,ijm1,2,1) * zv_ice(iim1,ijm1) 
    567                   ze22 = - akappa(iim1,ijm1,2,2) * zv_ice(iim1,ijm1) + akappa(iim1,ijm1,2,1) * zu_ice(iim1,ijm1) 
    568                   ze21 = - akappa(iim1,ijm1,1,1) * zv_ice(iim1,ijm1) - akappa(iim1,ijm1,1,2) * zu_ice(iim1,ijm1) 
    569                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    570                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    571                   zvis12 =       zviseta (iim1,ijm1) + dm 
    572                   zvis21 =       zviseta (iim1,ijm1) 
    573  
    574                   zdiag = zvis22 * ( ze11 + ze22 ) 
    575                   zs11(ji,jj,1,2) =  zs11(ji,jj,1,2) + zvis11 * ze11 + zdiag  
    576                   zs12(ji,jj,1,2) =  zs12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21 
    577                   zs22(ji,jj,1,2) =  zs22(ji,jj,1,2) + zvis11 * ze22 + zdiag 
    578                   zs21(ji,jj,1,2) =  zs21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12 
    579  
    580 #if defined key_agrif 
    581              END DO 
    582           END DO 
    583  
    584           DO jj = k_j1+1, k_jpj-1 
    585              DO ji = 2, jpim1 
    586 #endif 
    587                   zd1(ji,jj) =   & 
    588                      + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)  & 
    589                      - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)  & 
    590                      - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)  & 
    591                      + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)  & 
    592                      + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)  & 
    593                      + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)  & 
    594                      - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)  & 
    595                      - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 
    596                   zd2(ji,jj) =   & 
    597                      + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)  & 
    598                      - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)  & 
    599                      - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)  & 
    600                      + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)  & 
    601                      - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)  & 
    602                      - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)  & 
    603                      + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)  & 
    604                      + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 
    605                END DO 
     532            zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) ) 
     533!!!! this should be faster on scalar processor 
     534!           zresm = MAXVAL(  MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ),   & 
     535!              &                  ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) )  ) 
     536!!!! 
     537            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
     538 
     539            DO jj = k_j1, k_jpj 
     540               zu_a(:,jj) = zu_n(:,jj) 
     541               zv_a(:,jj) = zv_n(:,jj) 
    606542            END DO 
    607543 
    608             DO jj = k_j1+1, k_jpj-1 
    609                DO ji = 2, jpim1 
    610                   zunw = (  ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj)        & 
    611                      &    + ( za2(ji,jj) + zd2(ji,jj) ) * zc1(ji,jj) ) * zden(ji,jj) 
    612  
    613                   zvnw = (  ( za2(ji,jj) + zd2(ji,jj) ) * zb1(ji,jj)        & 
    614                      &    - ( za1(ji,jj) + zd1(ji,jj) ) * zb2(ji,jj) ) * zden(ji,jj) 
    615  
    616                   zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 
    617  
    618                   u_ice(ji,jj) = ( u_ice(ji,jj) + om * ( zunw - u_ice(ji,jj) ) * tmu(ji,jj) ) * zmask 
    619                   v_ice(ji,jj) = ( v_ice(ji,jj) + om * ( zvnw - v_ice(ji,jj) ) * tmu(ji,jj) ) * zmask 
    620                END DO 
     544            IF( zresm <= resl )   EXIT   iflag 
     545 
     546            !                                                   ! ================ ! 
     547         END DO    iflag                                        !  end Relaxation  ! 
     548         !                                                      ! ================ ! 
     549 
     550         IF( zindu == 0 ) THEN      ! even iteration 
     551            DO jj = k_j1 , k_jpj-1 
     552               zu0(:,jj) = zu_a(:,jj) 
     553               zv0(:,jj) = zv_a(:,jj) 
    621554            END DO 
    622  
    623             CALL lbc_lnk( u_ice, 'I', -1. ) 
    624             CALL lbc_lnk( v_ice, 'I', -1. ) 
    625  
    626             !---  5.2.5.4. Convergence test. 
    627             DO jj = k_j1+1 , k_jpj-1 
    628                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    629             END DO 
    630             zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 
    631             IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
    632  
    633             IF ( zresm <= resl) EXIT iflag 
    634  
    635          END DO iflag 
    636  
    637          zindu1 = 1.0 - zindu 
    638          DO jj = k_j1 , k_jpj-1 
    639             zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj) 
    640             zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj) 
    641          END DO 
    642       !                                                   ! ==================== ! 
     555         ENDIF 
     556         !                                                ! ==================== ! 
    643557      END DO                                              !  end loop over iter  ! 
    644558      !                                                   ! ==================== ! 
     559 
     560      ui_ice(:,:) = zu_a(:,1:jpj) 
     561      vi_ice(:,:) = zv_a(:,1:jpj) 
    645562 
    646563      IF(ln_ctl) THEN 
    647564         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter 
    648565         CALL prt_ctl_info(charout) 
    649          CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 
     566         CALL prt_ctl(tab2d_1=ui_ice, clinfo1=' lim_rhg  : ui_ice :', tab2d_2=vi_ice, clinfo2=' vi_ice :') 
    650567      ENDIF 
    651568 
  • trunk/NEMO/LIM_SRC_2/limrst_2.F90

    r823 r888  
    1717   !!---------------------------------------------------------------------- 
    1818   USE ice_2 
    19    USE dom_oce 
    20    USE ice_oce         ! ice variables 
     19   USE sbc_oce 
     20   USE sbc_ice 
    2121   USE daymod 
    2222 
     
    2727   PRIVATE 
    2828 
    29    PUBLIC   lim_rst_opn_2     ! routine called by ??? module 
    30    PUBLIC   lim_rst_write_2   ! routine called by ??? module 
    31    PUBLIC   lim_rst_read_2    ! routine called by ??? module 
     29   PUBLIC   lim_rst_opn_2     ! routine called by sbcice_lim_2.F90 
     30   PUBLIC   lim_rst_write_2   ! routine called by sbcice_lim_2.F90 
     31   PUBLIC   lim_rst_read_2    ! routine called by iceini_2.F90 
    3232 
    3333   LOGICAL, PUBLIC ::   lrst_ice         !: logical to control the ice restart write  
     
    3636   !!---------------------------------------------------------------------- 
    3737   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
    38    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limrst.F90,v 1.15 2007/06/29 14:54:06 opalod Exp $  
     38   !! $ Id: $ 
    3939   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4040   !!---------------------------------------------------------------------- 
     
    5757       
    5858      ! to get better performances with NetCDF format: 
    59       ! we open and define the ice restart file one ice time step before writing the data (-> at nitrst - 2*nfice + 1) 
    60       ! except if we write ice restart files every ice time step or if an ice restart file was writen at nitend - 2*nfice + 1 
    61       IF( kt == nitrst - 2*nfice + 1 .OR. nstock == nfice .OR. ( kt == nitend - nfice + 1 .AND. .NOT. lrst_ice ) ) THEN 
     59      ! we open and define the ice restart file one ice time step before writing the data (-> at nitrst - 2*nn_fsbc + 1) 
     60      ! except if we write ice restart files every ice time step or if an ice restart file was writen at nitend - 2*nn_fsbc + 1 
     61      IF( kt == nitrst - 2*nn_fsbc + 1 .OR. nstock == nn_fsbc .OR. ( kt == nitend - nn_fsbc + 1 .AND. .NOT. lrst_ice ) ) THEN 
    6262         ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    6363         IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     
    7272            CASE DEFAULT         ;   WRITE(numout,*) '             open ice restart NetCDF file: '//clname 
    7373            END SELECT 
    74             IF( kt == nitrst - 2*nfice + 1 ) THEN    
    75                WRITE(numout,*)         '             kt = nitrst - 2*nfice + 1 = ', kt,' date= ', ndastp 
    76             ELSE   ;   WRITE(numout,*) '             kt = '                       , kt,' date= ', ndastp 
     74            IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN    
     75               WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp 
     76            ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp 
    7777            ENDIF 
    7878         ENDIF 
     
    9090      !! ** purpose  :   output of sea-ice variable in a netcdf file 
    9191      !!---------------------------------------------------------------------- 
    92       INTEGER, INTENT(in) ::   kt     ! number of iteration 
    93       !! 
    94       INTEGER ::   iter     ! kt + nfice -1 
    95       !!---------------------------------------------------------------------- 
    96  
    97       iter = kt + nfice - 1   ! ice restarts are written at kt == nitrst - nfice + 1 
     92      INTEGER, INTENT(in) ::   kt   ! number of iteration 
     93      ! 
     94      INTEGER ::   iter   ! kt + nn_fsbc -1 
     95      !!---------------------------------------------------------------------- 
     96 
     97      iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1 
    9898 
    9999      IF( iter == nitrst ) THEN 
    100100         IF(lwp) WRITE(numout,*) 
    101101         IF(lwp) WRITE(numout,*) 'lim_rst_write_2 : write ice restart file  kt =', kt 
    102          IF(lwp) WRITE(numout,*) '~~~~~~~'          
     102         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    103103      ENDIF 
    104104 
     
    106106      ! ------------------  
    107107      !                                                                     ! calendar control 
    108       CALL iom_rstput( iter, nitrst, numriw, 'nfice' , REAL( nfice, wp) )      ! time-step  
    109       CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter , wp) )      ! date 
     108      CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) )  
    110109       
    111110      CALL iom_rstput( iter, nitrst, numriw, 'hicif' , hicif (:,:)   )      ! prognostic variables  
     
    119118      CALL iom_rstput( iter, nitrst, numriw, 'tbif2' , tbif  (:,:,2) ) 
    120119      CALL iom_rstput( iter, nitrst, numriw, 'tbif3' , tbif  (:,:,3) ) 
    121       CALL iom_rstput( iter, nitrst, numriw, 'u_ice' , u_ice (:,:)   ) 
    122       CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice (:,:)   ) 
    123       CALL iom_rstput( iter, nitrst, numriw, 'gtaux' , gtaux (:,:)   ) 
    124       CALL iom_rstput( iter, nitrst, numriw, 'gtauy' , gtauy (:,:)   ) 
     120      CALL iom_rstput( iter, nitrst, numriw, 'ui_ice', ui_ice(:,:)   ) 
     121      CALL iom_rstput( iter, nitrst, numriw, 'vi_ice', vi_ice(:,:)   ) 
    125122      CALL iom_rstput( iter, nitrst, numriw, 'qstoif', qstoif(:,:)   ) 
    126123      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq (:,:)   ) 
     
    175172      !! ** purpose  :   read of sea-ice variable restart in a netcdf file 
    176173      !!---------------------------------------------------------------------- 
    177       ! 
    178       REAL(wp) ::   zfice, ziter 
     174      REAL(wp) ::   ziter 
    179175      !!---------------------------------------------------------------------- 
    180176 
     
    182178         WRITE(numout,*) 
    183179         WRITE(numout,*) 'lim_rst_read_2 : read ice NetCDF restart file' 
    184          WRITE(numout,*) '~~~~~~~~' 
     180         WRITE(numout,*) '~~~~~~~~~~~~~~' 
    185181      ENDIF 
    186182 
    187183      CALL iom_open ( 'restart_ice_in', numrir, kiolib = jprstlib ) 
    188184 
    189       CALL iom_get( numrir, 'nfice' , zfice ) 
    190       CALL iom_get( numrir, 'kt_ice', ziter )     
    191       IF(lwp) WRITE(numout,*) '   read ice restart file at time step    : ', ziter 
     185      CALL iom_get( numrir, 'kt_ice' , ziter ) 
     186      IF(lwp) WRITE(numout,*) '   read ice restart file at time step    : ', INT( ziter ) 
    192187      IF(lwp) WRITE(numout,*) '   in any case we force it to nit000 - 1 : ', nit000 - 1 
    193188 
     
    196191      IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 )   & 
    197192         &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart',  & 
    198          &                   '   verify the file or rerun with the value 0 for the',        & 
    199          &                   '   control of time parameter  nrstdt' ) 
    200       IF( INT(zfice) /= nfice          .AND. ABS( nrstdt ) == 1 )   & 
    201          &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nfice in ice restart',  & 
    202193         &                   '   verify the file or rerun with the value 0 for the',        & 
    203194         &                   '   control of time parameter  nrstdt' ) 
     
    213204      CALL iom_get( numrir, jpdom_autoglo, 'tbif2' , tbif(:,:,2) )     
    214205      CALL iom_get( numrir, jpdom_autoglo, 'tbif3' , tbif(:,:,3) )     
    215       CALL iom_get( numrir, jpdom_autoglo, 'u_ice' , u_ice  )     
    216       CALL iom_get( numrir, jpdom_autoglo, 'v_ice' , v_ice  )     
    217       CALL iom_get( numrir, jpdom_autoglo, 'gtaux' , gtaux  )     
    218       CALL iom_get( numrir, jpdom_autoglo, 'gtauy' , gtauy  )     
     206      CALL iom_get( numrir, jpdom_autoglo, 'ui_ice', ui_ice )     
     207      CALL iom_get( numrir, jpdom_autoglo, 'vi_ice', vi_ice )     
    219208      CALL iom_get( numrir, jpdom_autoglo, 'qstoif', qstoif )     
    220209      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq  )     
  • trunk/NEMO/LIM_SRC_2/limtab_2.F90

    r823 r888  
    2121   !!---------------------------------------------------------------------- 
    2222   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    23    !! $Header$  
    24    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     23   !! $ Id: $ 
     24   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    2525   !!---------------------------------------------------------------------- 
    2626CONTAINS 
  • trunk/NEMO/LIM_SRC_2/limthd_2.F90

    r823 r888  
    44   !!              LIM thermo ice model : ice thermodynamic 
    55   !!====================================================================== 
     6   !! History :  1.0  !  00-01 (LIM) 
     7   !!            2.0  !  02-07 (C. Ethe, G. Madec) F90 
     8   !!            2.0  !  03-08 (C. Ethe)  add lim_thd_init 
     9   !!--------------------------------------------------------------------- 
    610#if defined key_lim2 
    711   !!---------------------------------------------------------------------- 
     
    1822   USE ice_2           ! LIM sea-ice variables 
    1923   USE ice_oce         ! sea-ice/ocean variables 
    20    USE flx_oce         ! sea-ice/ocean forcings variables  
     24   USE sbc_oce         !  
     25   USE sbc_ice         !  
    2126   USE thd_ice_2       ! LIM thermodynamic sea-ice variables 
    2227   USE dom_ice_2       ! LIM sea-ice domain 
     
    3035   PRIVATE 
    3136 
    32    !! * Routine accessibility 
    33    PUBLIC lim_thd_2       ! called by lim_step_2 
    34  
    35    !! * Module variables 
    36    REAL(wp)  ::            &  ! constant values 
    37       epsi20 = 1.e-20   ,  & 
    38       epsi16 = 1.e-16   ,  & 
    39       epsi04 = 1.e-04   ,  & 
    40       zzero  = 0.e0     ,  & 
    41       zone   = 1.e0 
     37   PUBLIC   lim_thd_2  ! called by lim_step 
     38 
     39   REAL(wp)  ::   epsi20 = 1.e-20   ,  &  ! constant values 
     40      &           epsi16 = 1.e-16   ,  & 
     41      &           epsi04 = 1.e-04   ,  & 
     42      &           zzero  = 0.e0     ,  & 
     43      &           zone   = 1.e0 
    4244 
    4345   !! * Substitutions 
     
    4648   !!-------- ------------------------------------------------------------- 
    4749   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    48    !! $Header$  
    49    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     50   !! $ Id: $ 
     51   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5052   !!---------------------------------------------------------------------- 
    5153 
     
    6870      !!             - back to the geographic grid 
    6971      !! 
    70       !! ** References : 
    71       !!       H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 
    72       !! 
    73       !! History : 
    74       !!   1.0  !  00-01 (LIM) 
    75       !!   2.0  !  02-07 (C. Ethe, G. Madec) F90 
     72      !! References :   Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 
    7673      !!--------------------------------------------------------------------- 
    7774      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    78  
     75      !! 
    7976      INTEGER  ::   ji, jj,    &   ! dummy loop indices 
    8077         nbpb  ,               &   ! nb of icy pts for thermo. cal. 
     
    9289         zfontn             ,  &   ! heat flux from snow thickness 
    9390         zfntlat, zpareff          ! test. the val. of lead heat budget 
    94       REAL(wp), DIMENSION(jpi,jpj) :: & 
    95          zhicifp            ,  &   ! ice thickness for outputs 
    96          zqlbsbq                   ! link with lead energy budget qldif 
    97       REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 
    98          zmsk                      ! working array 
     91      REAL(wp), DIMENSION(jpi,jpj) ::   zhicifp,   &  ! ice thickness for outputs 
     92         &                              zqlbsbq       ! link with lead energy budget qldif 
     93      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmsk      ! working array 
    9994      !!------------------------------------------------------------------- 
    10095 
    101       IF( kt == nit000  )   CALL lim_thd_init_2  ! Initialization (first time-step only) 
     96      IF( kt == nit000 )   CALL lim_thd_init_2  ! Initialization (first time-step only) 
    10297    
    10398      !-------------------------------------------! 
     
    173168      !-------------------------------------------------------------------------- 
    174169 
     170      sst_m(:,:) = sst_m(:,:) + rt0 
     171 
    175172      !CDIR NOVERRCHK 
    176173      DO jj = 1, jpj 
     
    188185            !  temperature and turbulent mixing (McPhee, 1992) 
    189186            zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity 
    190             fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_io(ji,jj) - tfu(ji,jj) )  
     187            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) )  
    191188            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice 
    192189                         
    193190            !  partial computation of the lead energy budget (qldif) 
    194191            zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting 
    195             zfnsol         = qnsr_oce(ji,jj)  !  total non solar flux 
    196             qldif(ji,jj)   = tms(ji,jj) * ( qsr_oce(ji,jj) * ( 1.0 - thcm(ji,jj) )   & 
     192            zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean 
     193            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   & 
    197194               &                               + zfnsol + fdtcn(ji,jj) - zfontn     & 
    198195               &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   & 
     
    206203             
    207204            !  energy needed to bring ocean surface layer until its freezing 
    208             qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_io(ji,jj) ) * ( 1 - zinda ) 
     205            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda ) 
    209206             
    210207            !  calculate oceanic heat flux. 
     
    216213      END DO 
    217214       
     215      sst_m(:,:) = sst_m(:,:) - rt0 
    218216       
    219217      !         Select icy points and fulfill arrays for the vectorial grid. 
     
    258256         CALL tab_2d_1d_2( nbpb, fr1_i0_1d  (1:nbpb)     , fr1_i0     , jpi, jpj, npb(1:nbpb) ) 
    259257         CALL tab_2d_1d_2( nbpb, fr2_i0_1d  (1:nbpb)     , fr2_i0     , jpi, jpj, npb(1:nbpb) ) 
    260          CALL tab_2d_1d_2( nbpb, qnsr_ice_1d(1:nbpb)     , qnsr_ice   , jpi, jpj, npb(1:nbpb) ) 
     258         CALL tab_2d_1d_2( nbpb, qns_ice_1d (1:nbpb)     , qns_ice    , jpi, jpj, npb(1:nbpb) ) 
    261259#if ! defined key_coupled 
    262260         CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb)     , qla_ice    , jpi, jpj, npb(1:nbpb) ) 
     
    404402         CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif  : ', tab2d_2=fsbbq , clinfo2=' fsbbq  : ') 
    405403      ENDIF 
    406  
     404       ! 
    407405    END SUBROUTINE lim_thd_2 
    408406 
     
    419417      !! 
    420418      !! ** input   :   Namelist namicether 
    421       !! 
    422       !! history : 
    423       !!  8.5  ! 03-08 (C. Ethe) original code 
    424419      !!------------------------------------------------------------------- 
    425420      NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax  ,        & 
  • trunk/NEMO/LIM_SRC_2/limthd_lac_2.F90

    r823 r888  
    3030   !!---------------------------------------------------------------------- 
    3131   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    32    !! $Header$  
     32   !! $ Id: $ 
    3333   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3434   !!---------------------------------------------------------------------- 
  • trunk/NEMO/LIM_SRC_2/limthd_zdf_2.F90

    r823 r888  
    44   !!                thermodynamic growth and decay of the ice  
    55   !!====================================================================== 
     6   !! History :  1.0  !  01-04 (LIM) Original code 
     7   !!            2.0  !  02-08 (C. Ethe, G. Madec) F90 
     8   !!---------------------------------------------------------------------- 
    69#if defined key_lim2 
    710   !!---------------------------------------------------------------------- 
    811   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
     12   !!---------------------------------------------------------------------- 
    913   !!---------------------------------------------------------------------- 
    1014   !!   lim_thd_zdf_2 : vertical accr./abl. and lateral ablation of sea ice 
     
    2226   PRIVATE 
    2327 
    24    !! * Routine accessibility 
    25    PUBLIC lim_thd_zdf_2      ! called by lim_thd_2 
    26  
    27    !! * Module variables 
    28    REAL(wp)  ::           &  ! constant values 
    29       epsi20 = 1.e-20  ,  & 
    30       epsi13 = 1.e-13  ,  & 
    31       zzero  = 0.e0    ,  & 
    32       zone   = 1.e0 
     28   PUBLIC   lim_thd_zdf_2        ! called by lim_thd_2 
     29 
     30   REAL(wp) ::   epsi20 = 1.e-20  ,  &  ! constant values 
     31      &          epsi13 = 1.e-13  ,  & 
     32      &          zzero  = 0.e0    ,  & 
     33      &          zone   = 1.e0 
    3334   !!---------------------------------------------------------------------- 
    3435   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    35    !! $Header$  
    36    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    37    !!---------------------------------------------------------------------- 
     36   !! $ Id: $ 
     37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     38   !!---------------------------------------------------------------------- 
     39 
    3840CONTAINS 
    3941 
     
    6466      !!              - Performs lateral ablation 
    6567      !! 
    66       !! References : 
    67       !!   Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646    
    68       !!   Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268   
     68      !! References : Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646    
     69      !!              Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268   
     70      !!------------------------------------------------------------------ 
     71      INTEGER, INTENT(in) ::   kideb    ! Start point on which the  the computation is applied 
     72      INTEGER, INTENT(in) ::   kiut     ! End point on which the  the computation is applied 
    6973      !! 
    70       !! History : 
    71       !!   original    : 01-04 (LIM) 
    72       !!   addition    : 02-08 (C. Ethe, G. Madec) 
    73       !!------------------------------------------------------------------ 
    74       !! * Arguments 
    75       INTEGER , INTENT (in) ::  & 
    76          kideb ,  &  ! Start point on which the  the computation is applied 
    77          kiut        ! End point on which the  the computation is applied 
    78  
    79       !! * Local variables 
    8074      INTEGER ::   ji       ! dummy loop indices 
    81  
    82       REAL(wp) , DIMENSION(jpij,2) ::  & 
    83          zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 ) 
    84  
     75      REAL(wp), DIMENSION(jpij,2) ::   zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 ) 
    8576      REAL(wp), DIMENSION(jpij) ::  & 
    8677         ztsmlt      &    ! snow/ice surface melting temperature 
     
    9788         , zts_old   &    ! previous surface temperature 
    9889         , zidsn , z1midsn , zidsnic ! tempory variables 
    99  
    100       REAL(wp), DIMENSION(jpij) :: & 
     90      REAL(wp), DIMENSION(jpij) ::   & 
    10191          zfnet       &  ! net heat flux at the top surface( incl. conductive heat flux) 
    10292          , zsprecip  &    ! snow accumulation 
     
    10999          , zfrld_1d    &    ! new sea/ice fraction 
    110100          , zep            ! internal temperature of the 2nd layer of the snow/ice system 
    111  
    112101       REAL(wp), DIMENSION(3) :: &  
    113102          zplediag  &    ! principle diagonal, subdiag. and supdiag. of the  
     
    115104          , zsupdiag  &    ! of the temperatures inside the snow-ice system 
    116105          , zsmbr          ! second member 
    117  
    118106       REAL(wp) :: &  
    119107          zhsu     &     ! thickness of surface layer 
     
    130118          , zb2 , zd2 , zb3 , zd3 & 
    131119          , ztint          ! equivalent temperature at the snow-ice interface 
    132  
    133120       REAL(wp) :: &  
    134121          zexp      &     ! exponential function of the ice thickness 
     
    148135          , zdtic        &  ! ice internal temp. increment 
    149136          , zqnes          ! conductive energy due to ice melting in the first ice layer 
    150  
    151137       REAL(wp) :: &  
    152138          ztbot     &      ! temperature at the bottom surface 
     
    162148          , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables 
    163149          , ztb2, ztb3 
    164  
    165150       REAL(wp) :: &  
    166151          zdrmh         &   ! change in snow/ice thick. after snow-ice formation 
     
    181166       !     Computation of energies due to surface and bottom melting  
    182167       !----------------------------------------------------------------------- 
    183         
    184168        
    185169       DO ji = kideb , kiut 
     
    201185       END DO 
    202186 
    203  
    204187       !------------------------------------------- 
    205188       !  2. Calculate some intermediate variables.   
     
    265248       !     - qstbif_1d, total energy stored in brine pockets (updating) 
    266249       !------------------------------------------------------------------- 
    267  
    268250 
    269251       DO ji = kideb , kiut 
     
    288270       END DO 
    289271 
    290  
    291272       !-------------------------------------------------------------------------------- 
    292273       !  4. Computation of the surface temperature : determined by considering the  
     
    333314          !---computation of the energy balance function  
    334315          zfts    = - z1mi0 (ji) * qsr_ice_1d(ji)   & ! net absorbed solar radiation 
    335              &      - qnsr_ice_1d(ji)                & ! total non solar radiation 
    336              &      - zfcsu (ji)                  ! conductive heat flux from the surface 
     316             &      - qns_ice_1d(ji)                & ! total non solar radiation 
     317             &      - zfcsu (ji)                      ! conductive heat flux from the surface 
    337318          !---computation of surface temperature increment   
    338319          zdts    = -zfts / zdfts 
     
    360341          sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 
    361342#if ! defined key_coupled 
    362           qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
     343          qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
    363344          qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
    364345#endif 
     
    366347       END DO 
    367348 
    368       
    369  
    370349       !     5.2. Calculate available heat for surface ablation.  
    371350       !--------------------------------------------------------------------- 
    372351 
    373352       DO ji = kideb, kiut 
    374           zfnet(ji) = qnsr_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)           
     353          zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)           
    375354          zfnet(ji) = MAX( zzero , zfnet(ji) ) 
    376355          zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) ) 
     
    730709          dvnbq_1d(ji) = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) 
    731710          dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn 
    732           !  case of natural freshwater flux 
    733 #if defined key_lim_fdd   
    734           rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d(ji) )  * rhoic 
     711          !---  volume change of ice and snow (used for ocean-ice freshwater flux computation) 
     712          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic 
    735713          rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn 
    736714 
    737 #else 
    738           rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * (  ( zhicnew - h_ice_1d(ji) ) * rhoic   & 
    739              &                                                    + ( zhsnnew - h_snow_1d(ji) ) * rhosn ) 
    740 #endif 
    741  
    742715          !---  Actualize new snow and ice thickness. 
    743  
    744716          h_snow_1d(ji)  = zhsnnew 
    745           h_ice_1d(ji)  = zhicnew 
     717          h_ice_1d (ji)  = zhicnew 
    746718 
    747719       END DO 
     
    799771          qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji) 
    800772          frld_1d(ji)    = zfrld_1d(ji) 
    801  
    802        END DO 
    803         
     773          ! 
     774       END DO 
     775       !  
    804776    END SUBROUTINE lim_thd_zdf_2 
     777 
    805778#else 
    806    !!====================================================================== 
    807    !!                       ***  MODULE limthd_zdf_2   *** 
    808    !!                              no sea ice model   
    809    !!====================================================================== 
     779   !!---------------------------------------------------------------------- 
     780   !!   Default Option                                     NO sea-ice model 
     781   !!---------------------------------------------------------------------- 
    810782CONTAINS 
    811783   SUBROUTINE lim_thd_zdf_2          ! Empty routine 
    812784   END SUBROUTINE lim_thd_zdf_2 
    813785#endif 
    814  END MODULE limthd_zdf_2 
     786 
     787   !!====================================================================== 
     788END MODULE limthd_zdf_2 
  • trunk/NEMO/LIM_SRC_2/limtrp_2.F90

    r823 r888  
    3030 
    3131   !! * Routine accessibility 
    32    PUBLIC lim_trp_2     ! called by ice_step_2 
     32   PUBLIC lim_trp_2     ! called by sbc_ice_lim_2 
    3333 
    3434   !! * Shared module variables 
     
    4848   !!---------------------------------------------------------------------- 
    4949   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    50    !! $Header$  
     50   !! $ Id: $ 
    5151   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    5252   !!---------------------------------------------------------------------- 
     
    112112         DO jj = 1, jpjm1 
    113113            DO ji = 1, jpim1 
    114                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 ) ) 
    115                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 ) ) 
     114               zui_u(ji,jj) = ( ui_ice(ji+1,jj  ) + ui_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) ) 
     115               zvi_v(ji,jj) = ( vi_ice(ji  ,jj+1) + vi_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) ) 
    116116            END DO 
    117117         END DO 
     
    128128         IF (lk_mpp ) CALL mpp_max(zcfl) 
    129129 
    130          IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl 
     130         IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl 
    131131 
    132132         ! content of properties 
  • trunk/NEMO/LIM_SRC_2/limwri_2.F90

    r823 r888  
    99#if defined key_lim2 
    1010   !!---------------------------------------------------------------------- 
    11    !!   'key_lim2'  i                                 LIM 2.0 sea-ice model 
     11   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
    1313   !!---------------------------------------------------------------------- 
     
    1515   !!   lim_wri_init_2 : initialization and namelist read 
    1616   !!---------------------------------------------------------------------- 
    17    USE ioipsl 
    18    USE dianam    ! build name of file (routine) 
    1917   USE phycst 
    2018   USE dom_oce 
    2119   USE daymod 
    22    USE in_out_manager 
    2320   USE ice_oce         ! ice variables 
    24    USE flx_oce 
     21   USE sbc_oce 
     22   USE sbc_ice 
    2523   USE dom_ice_2 
    2624   USE ice_2 
     25 
    2726   USE lbclnk 
     27   USE dianam          ! build name of file (routine) 
     28   USE in_out_manager 
     29   USE ioipsl 
    2830 
    2931   IMPLICIT NONE 
    3032   PRIVATE 
    3133 
    32    PUBLIC   lim_wri_2        ! routine called by lim_step_2.F90 
     34   PUBLIC   lim_wri_2      ! routine called by sbc_ice_lim_2 
    3335 
    3436   INTEGER, PARAMETER                       ::   jpnoumax = 40   ! maximum number of variable for ice output 
     
    4951      zone   = 1.e0 
    5052 
    51    !!---------------------------------------------------------------------- 
    52    !!  LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    53    !! $Header$ 
     53   !! * Substitutions 
     54#   include "vectopt_loop_substitute.h90" 
     55   !!---------------------------------------------------------------------- 
     56   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006) 
     57   !! $ Id: $ 
    5458   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5559   !!---------------------------------------------------------------------- 
     
    7983      !!------------------------------------------------------------------- 
    8084      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    81  
     85      !! 
    8286      INTEGER  ::   ji, jj, jf                      ! dummy loop indices 
    8387      CHARACTER(len = 40)  ::   clhstnam, clop 
     
    9094 
    9195      !                                          !--------------------! 
    92       IF ( kt == nit000 ) THEN                !   Initialisation   ! 
     96      IF( kt == nit000 ) THEN                    !   Initialisation   ! 
    9397         !                                       !--------------------! 
    9498         CALL lim_wri_init_2  
     
    97101!!Chris         clop     = "ave(only(x))"      !ibug  namelist parameter a ajouter 
    98102         clop     = "ave(x)" 
    99          zout     = nwrite * rdt_ice / nfice 
     103         zout     = nwrite * rdt_ice / nn_fsbc 
    100104         zsec     = 0. 
    101105         niter    = 0 
     
    110114          
    111115         DO jf = 1, noumef 
    112             IF ( nc(jf) == 1 )   CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj   & 
    113                   &                                , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout ) 
     116            IF( nc(jf) == 1 )   CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj   & 
     117                                               , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    114118         END DO 
    115119         CALL histend( nice ) 
    116           
     120         ! 
    117121      ENDIF 
    118122      !                                          !--------------------! 
     
    120124      !                                          !--------------------! 
    121125 
    122 !!gm  change the print below to have it only at output time step or when nitend =< 100 
    123       IF(lwp) THEN 
    124          WRITE(numout,*) 
    125          WRITE(numout,*) 'lim_wri_2 : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, kt + nfice - 1 
    126          WRITE(numout,*) '~~~~~~~~~ ' 
    127       ENDIF 
    128  
    129       !-- calculs des valeurs instantanees 
     126      !-- Store instantaneous values in zcmo 
    130127       
    131128      zcmo(:,:, 1:jpnoumax ) = 0.e0  
    132129      DO jj = 2 , jpjm1 
    133          DO ji = 2 , jpim1 
     130         DO ji = fs_2 , fs_jpim1 
    134131            zindh  = MAX( zzero , SIGN( zone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) ) 
    135132            zinda  = MAX( zzero , SIGN( zone , ( 1.0 - frld(ji,jj) ) - 0.10 ) ) 
     
    142139            zcmo(ji,jj,5)  = sist  (ji,jj) 
    143140            zcmo(ji,jj,6)  = fbif  (ji,jj) 
    144             zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
    145                                       + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
     141            zcmo(ji,jj,7)  = zindb * (  ui_ice(ji,jj  ) * tmu(ji,jj  ) + ui_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
     142                                      + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
    146143                                  / ztmu  
    147144 
    148             zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
    149                                       + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
     145            zcmo(ji,jj,8)  = zindb * (  vi_ice(ji,jj  ) * tmu(ji,jj  ) + vi_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
     146                                      + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
    150147                                  / ztmu 
    151             zcmo(ji,jj,9)  = sst_io(ji,jj) 
    152             zcmo(ji,jj,10) = sss_io(ji,jj) 
    153  
    154             zcmo(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj) 
    155             zcmo(ji,jj,12) = fsolar (ji,jj) 
    156             zcmo(ji,jj,13) = fnsolar(ji,jj) 
     148            zcmo(ji,jj,9)  = sst_m(ji,jj) 
     149            zcmo(ji,jj,10) = sss_m(ji,jj) 
     150            zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj) 
     151            zcmo(ji,jj,12) = qsr(ji,jj) 
     152            zcmo(ji,jj,13) = qns(ji,jj) 
    157153            ! See thersf for the coefficient 
    158             zcmo(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_io(ji,jj) + epsi16 ) / soce 
    159             zcmo(ji,jj,15) = gtaux(ji,jj) 
    160             zcmo(ji,jj,16) = gtauy(ji,jj) 
    161             zcmo(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice (ji,jj) + frld(ji,jj) * qsr_oce (ji,jj) 
    162             zcmo(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qnsr_ice(ji,jj) + frld(ji,jj) * qnsr_oce(ji,jj) 
     154            zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ??? 
     155            zcmo(ji,jj,15) = utaui_ice(ji,jj) 
     156            zcmo(ji,jj,16) = vtaui_ice(ji,jj) 
     157            zcmo(ji,jj,17) = qsr_ice(ji,jj) 
     158            zcmo(ji,jj,18) = qns_ice(ji,jj) 
    163159            zcmo(ji,jj,19) = sprecip(ji,jj) 
    164160         END DO 
     
    175171         END DO 
    176172          
    177          IF ( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN  
     173         IF( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN 
    178174            CALL lbc_lnk( zfield, 'T', -1. ) 
    179175         ELSE  
     
    181177         ENDIF 
    182178          
    183          IF ( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 
     179         IF( nc(jf) == 1 )  CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 
    184180          
    185181      END DO 
    186182       
    187       IF ( ( nfice * niter + nit000 - 1 ) >= nitend ) THEN 
    188          CALL histclo( nice )  
    189       ENDIF 
     183      IF( ( nn_fsbc * niter + nit000 - 1 ) >= nitend )   CALL histclo( nice )  
    190184      ! 
    191185   END SUBROUTINE lim_wri_2 
     
    225219         field_13, field_14, field_15, field_16, field_17, field_18,   & 
    226220         field_19 
    227 !!gm      NAMELIST/namiceout/ noumef, & 
    228 !!           zfield( 1), zfield( 2), zfield( 3), zfield( 4), zfield( 5),   & 
    229 !!           zfield( 6), zfield( 7), zfield( 8), zfield( 9), zfield(10),   & 
    230 !!           zfield(11), zfield(12), zfield(13), zfield(14), zfield(15),   & 
    231 !!gm         zfield(16), zfield(17), zfield(18), zfield(19) 
    232       !!------------------------------------------------------------------- 
    233  
    234       ! Read Namelist namicewri 
    235       REWIND ( numnam_ice ) 
     221      !!------------------------------------------------------------------- 
     222 
     223      REWIND ( numnam_ice )                ! Read Namelist namicewri 
    236224      READ   ( numnam_ice  , namiceout ) 
    237225       
    238       zfield(1) = field_1 
    239       zfield(2) = field_2 
    240       zfield(3) = field_3 
    241       zfield(4) = field_4 
    242       zfield(5) = field_5 
    243       zfield(6) = field_6 
    244       zfield(7) = field_7 
    245       zfield(8) = field_8 
    246       zfield(9) = field_9 
     226      zfield( 1) = field_1 
     227      zfield( 2) = field_2 
     228      zfield( 3) = field_3 
     229      zfield( 4) = field_4 
     230      zfield( 5) = field_5 
     231      zfield( 6) = field_6 
     232      zfield( 7) = field_7 
     233      zfield( 8) = field_8 
     234      zfield( 9) = field_9 
    247235      zfield(10) = field_10 
    248236      zfield(11) = field_11 
     
    274262         DO nf = 1 , noumef          
    275263            WRITE(numout,*) '   ', titn(nf), '   ', nam(nf),'      ', uni(nf),'  ', nc(nf),'        ', cmulti(nf),   & 
    276                '        ', cadd(nf) 
     264               &       '        ', cadd(nf) 
    277265         END DO 
    278266      ENDIF 
  • trunk/NEMO/LIM_SRC_2/limwri_dimg_2.h90

    r823 r888  
    22   !!---------------------------------------------------------------------- 
    33   !!  LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    4    !! $Header$ 
     4   !! $ Id: $ 
    55   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    66   !!---------------------------------------------------------------------- 
     
    8282 
    8383       zsto     = rdt_ice 
    84        zout     = nwrite * rdt_ice / nfice 
     84       zout     = nwrite * rdt_ice / nn_fsbc 
    8585       zsec     = 0. 
    8686       niter    = 0 
     
    106106          zcmo(ji,jj,5)  = sist  (ji,jj) 
    107107          zcmo(ji,jj,6)  = fbif  (ji,jj) 
    108           zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
    109                + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
     108          zcmo(ji,jj,7)  = zindb * (  ui_ice(ji,jj  ) * tmu(ji,jj  ) + ui_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
     109             &                      + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
    110110               / ztmu  
    111111 
    112           zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
    113                + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
     112          zcmo(ji,jj,8)  = zindb * (  vi_ice(ji,jj  ) * tmu(ji,jj  ) + vi_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
     113             &                      + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
    114114               / ztmu 
    115           zcmo(ji,jj,9)  = sst_io(ji,jj) 
    116           zcmo(ji,jj,10) = sss_io(ji,jj) 
    117  
    118           zcmo(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj) 
    119           zcmo(ji,jj,12) = fsolar (ji,jj) 
    120           zcmo(ji,jj,13) = fnsolar(ji,jj) 
     115          zcmo(ji,jj,9)  = sst_m(ji,jj) 
     116          zcmo(ji,jj,10) = sss_m(ji,jj) 
     117 
     118          zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj) 
     119          zcmo(ji,jj,12) = qsr(ji,jj) 
     120          zcmo(ji,jj,13) = qns(ji,jj) 
    121121          ! See thersf for the coefficient 
    122           zcmo(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_io(ji,jj) + epsi16 ) / soce 
    123           zcmo(ji,jj,15) = gtaux(ji,jj) 
    124           zcmo(ji,jj,16) = gtauy(ji,jj) 
    125           zcmo(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice (ji,jj) + frld(ji,jj) * qsr_oce (ji,jj) 
    126           zcmo(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qnsr_ice(ji,jj) + frld(ji,jj) * qnsr_oce(ji,jj) 
     122          zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     123          zcmo(ji,jj,15) = utaui_ice(ji,jj) 
     124          zcmo(ji,jj,16) = vtaui_ice(ji,jj) 
     125          zcmo(ji,jj,17) = qsr_ice(ji,jj) 
     126          zcmo(ji,jj,18) = qns_ice(ji,jj) 
    127127          zcmo(ji,jj,19) = sprecip(ji,jj) 
    128128       END DO 
     
    132132    nmoyice = nmoyice + 1  
    133133    ! compute mean value if it is time to write on file 
    134     IF ( MOD(kt+nfice-1-nit000+1,nwrite) == 0 ) THEN 
     134    IF ( MOD(kt+nn_fsbc-1-nit000+1,nwrite) == 0 ) THEN 
    135135       rcmoy(:,:,:) = rcmoy(:,:,:) / FLOAT(nmoyice) 
    136136#else   
    137        IF ( MOD(kt-nfice-1-nit000+1,nwrite) == 0 ) THEN  
     137       IF ( MOD(kt-nn_fsbc-1-nit000+1,nwrite) == 0 ) THEN  
    138138          !  case of instantaneaous output rcmoy(:,:, 1:jpnoumax ) = 0.e0 
    139139          DO jj = 2 , jpjm1 
     
    149149                rcmoy(ji,jj,5)  = sist  (ji,jj) 
    150150                rcmoy(ji,jj,6)  = fbif  (ji,jj) 
    151                 rcmoy(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
    152                      + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
     151                rcmoy(ji,jj,7)  = zindb * (  ui_ice(ji,jj  ) * tmu(ji,jj  ) + ui_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
     152                   &                       + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
    153153                     / ztmu 
    154154 
    155                 rcmoy(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
    156                      + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
     155                rcmoy(ji,jj,8)  = zindb * (  vi_ice(ji,jj  ) * tmu(ji,jj  ) + vi_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
     156                   &                       + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
    157157                     / ztmu 
    158                 rcmoy(ji,jj,9)  = sst_io(ji,jj) 
    159                 rcmoy(ji,jj,10) = sss_io(ji,jj) 
    160  
    161                 rcmoy(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj) 
    162                 rcmoy(ji,jj,12) = fsolar (ji,jj) 
    163                 rcmoy(ji,jj,13) = fnsolar(ji,jj) 
     158                rcmoy(ji,jj,9)  = sst_m(ji,jj) 
     159                rcmoy(ji,jj,10) = sss_m(ji,jj) 
     160 
     161                rcmoy(ji,jj,11) = qns(ji,jj) + qsr(ji,jj) 
     162                rcmoy(ji,jj,12) = qsr(ji,jj) 
     163                rcmoy(ji,jj,13) = qns(ji,jj) 
    164164                ! See thersf for the coefficient 
    165                 rcmoy(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_io(ji,jj) + epsi16 ) / soce 
    166                 rcmoy(ji,jj,15) = gtaux(ji,jj) 
    167                 rcmoy(ji,jj,16) = gtauy(ji,jj) 
    168                 rcmoy(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice (ji,jj) + frld(ji,jj) * qsr_oce (ji,jj) 
    169                 rcmoy(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qnsr_ice(ji,jj) + frld(ji,jj) * qnsr_oce(ji,jj) 
     165                rcmoy(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     166                rcmoy(ji,jj,15) = utaui_ice(ji,jj) 
     167                rcmoy(ji,jj,16) = vtaui_ice(ji,jj) 
     168                rcmoy(ji,jj,17) = qsr_ice(ji,jj) 
     169                rcmoy(ji,jj,18) = qns_ice(ji,jj) 
    170170                rcmoy(ji,jj,19) = sprecip(ji,jj) 
    171171             END DO 
     
    201201          rcmoy(:,:,:) = 0.0 
    202202          nmoyice = 0  
    203        END IF     !  MOD(kt+nfice-1-nit000+1, nwrite == 0 ) ! 
     203       END IF     !  MOD(kt+nn_fsbc-1-nit000+1, nwrite == 0 ) ! 
    204204 
    205205     END SUBROUTINE lim_wri_2 
  • trunk/NEMO/LIM_SRC_2/par_ice_2.F90

    r823 r888  
    77   !!---------------------------------------------------------------------- 
    88   !!  LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    9    !! $Header$ 
     9   !! $ Id: $ 
    1010   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    1111   !!---------------------------------------------------------------------- 
  • trunk/NEMO/LIM_SRC_2/thd_ice_2.F90

    r823 r888  
    99   !!---------------------------------------------------------------------- 
    1010   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    11    !! $Header$ 
     11   !! $ Id: $ 
    1212   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    1313   !!---------------------------------------------------------------------- 
     
    5757      fr1_i0_1d   ,     &  !:    "                  "      fr1_i0 
    5858      fr2_i0_1d   ,     &  !:    "                  "      fr2_i0 
    59       qnsr_ice_1d ,     &  !:    "                  "      qns_ice 
     59      qns_ice_1d ,     &  !:    "                  "      qns_ice 
    6060      qfvbq_1d    ,     &  !:    "                  "      qfvbq 
    6161      sist_1d     ,     &  !:    "                  "      sist 
Note: See TracChangeset for help on using the changeset viewer.