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

Changeset 717


Ignore:
Timestamp:
2007-10-16T13:03:55+02:00 (17 years ago)
Author:
smasson
Message:

finalize the first set of modifications related to ticket:3

Location:
trunk/NEMO
Files:
29 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMO/C1D_SRC/step1d.F90

    r714 r717  
    154154         CALL prt_ctl(tab2d_1=qns    , clinfo1=' qns  -   : ', mask1=tmask, ovlap=1) 
    155155         CALL prt_ctl(tab2d_1=qsr    , clinfo1=' qsr  -   : ', mask1=tmask, ovlap=1) 
    156          CALL prt_ctl(tab2d_1=runoff , clinfo1=' runoff   : ', mask1=tmask, ovlap=1) 
    157156         CALL prt_ctl(tab3d_1=tmask  , clinfo1=' tmask    : ', mask1=tmask, ovlap=1, kdim=jpk) 
    158157         CALL prt_ctl(tab3d_1=tn     , clinfo1=' sst  -   : ', mask1=tmask, ovlap=1, kdim=1) 

  • trunk/NEMO/LIM_SRC/ice.F90

    r699 r717  
    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_ice_lim 
    79   !!---------------------------------------------------------------------- 
    810   !!   'key_ice_lim' :                                   LIM 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    !! $Id$ 
    15    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    16    !!---------------------------------------------------------------------- 
    17    !! * Modules used 
    1812   USE par_ice          ! LIM sea-ice parameters 
    1913 
    … …  
    2115   PRIVATE 
    2216 
    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 
     17   !!* ice-dynamic namelist (namicedyn) * 
     18   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation 
     19   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation 
     20   REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic 
     21   REAL(wp), PUBLIC ::   alpha  = 0.5       !: coefficient for semi-implicit coriolis 
     22   REAL(wp), PUBLIC ::   dm     = 0.6e+03   !: diffusion constant for dynamics 
     23   REAL(wp), PUBLIC ::   om     = 0.5       !: relaxation constant 
     24   REAL(wp), PUBLIC ::   resl   = 5.0e-05   !: maximum value for the residual of relaxation 
     25   REAL(wp), PUBLIC ::   cw     = 5.0e-03   !: drag coefficient for oceanic stress 
     26   REAL(wp), PUBLIC ::   angvg  = 0.e0      !: turning angle for oceanic stress 
     27   REAL(wp), PUBLIC ::   pstar  = 1.0e+04   !: first bulk-rheology parameter 
     28   REAL(wp), PUBLIC ::   c_rhg  = 20.e0     !: second bulk-rhelogy parameter 
     29   REAL(wp), PUBLIC ::   etamn  = 0.e+07    !: minimun value for viscosity 
     30   REAL(wp), PUBLIC ::   creepl = 2.e-08    !: creep limit 
     31   REAL(wp), PUBLIC ::   ecc    = 2.e0      !: eccentricity of the elliptical yield curve 
     32   REAL(wp), PUBLIC ::   ahi0   = 350.e0    !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    2733 
    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) 
     34   REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc ) 
     35   REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw 
     36   REAL(wp), PUBLIC ::   sangvg, cangvg     !: sin and cos of the turning angle for ocean stress 
     37   REAL(wp), PUBLIC ::   pstarh             !: pstar / 2.0 
    4238 
    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 
     39   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
     40   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
     41   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses 
     42   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity 
    4843 
    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 
     44   !!* diagnostic quantities 
     45!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   firic         !: IR flux over the ice (only used for outputs) 
     46!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fcsic         !: Sensible heat flux over the ice (only used for outputs) 
     47!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fleic         !: Latent heat flux over the ice (only used for outputs) 
     48!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qlatic        !: latent flux (only used for outputs) 
     49   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvosif       !: Variation of volume at surface (only used for outputs) 
     50   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvobif       !: Variation of ice volume at the bottom ice (only used for outputs) 
     51   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdvolif       !: Total variation of ice volume (only used for outputs) 
     52   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvonif       !: Lateral Variation of ice volume (only used for outputs) 
    5553 
    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  
     54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin ??? degree ??? I don't know) 
     55   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS 
     56   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness 
     57   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnif         !: Snow thickness 
     58   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicifp        !: Ice production/melting 
     59   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   frld          !: Leads fraction = 1-a/totalarea 
     60   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   phicif        !: ice thickness  at previous time  
     61   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pfrld         !: Leads fraction at previous time   
     62   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qstoif        !: Energy stored in the brine pockets 
     63   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fbif          !: Heat flux at the ice base 
     64   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmsnif       !: Variation of snow mass 
     65   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass 
     66   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean) 
     67   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing  
     68   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean 
     69   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2) 
     70   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget 
     71   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice 
     72   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max heat contained in brine pockets (?) 
     73   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?) 
     74   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?) 
     75   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy in case of toral lateral ablation (?) 
     76   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice 
    5977 
    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 
     78   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   albege        !: Albedo of the snow or ice (only for outputs) 
     79   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   albecn        !: Albedo of the ocean (only for outputs) 
     80   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tauc          !: Cloud optical depth 
    9381 
    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) 
     82   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ui_ice, vi_ice   !: two components of the ice   velocity at I-point (m/s) 
     83   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ui_oce, vi_oce   !: two components of the ocean velocity at I-point (m/s) 
    9984 
     85   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpsmax)     ::   scal0   !: ??? 
     86   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer 
    10087 
    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) 
     88!! REAL(wp), DIMENSION(jpi,jpj,0:jpkmax+1) ::   reslum        !: Relative absorption of solar radiation in each ocean level 
    10489 
    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           !: 
     90   !!* moment used in the advection scheme 
     91   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume 
     92   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume 
     93   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxa,   sya,   sxxa,   syya,   sxya     !: for ice cover area 
     94   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc0,  syc0,  sxxc0,  syyc0,  sxyc0    !: for heat content of snow 
     95   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc1,  syc1,  sxxc1,  syyc1,  sxyc1    !: for heat content of 1st ice layer 
     96   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc2,  syc2,  sxxc2,  syyc2,  sxyc2    !: for heat content of 2nd ice layer 
     97   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxst,  syst,  sxxst,  syyst,  sxyst    !: for heat content of brine pockets 
    12298 
    12399#else 
    … …  
    127103#endif 
    128104 
     105   !!---------------------------------------------------------------------- 
     106   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006) 
     107   !! $Id$ 
     108   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    129109   !!====================================================================== 
    130110END MODULE ice 

  • trunk/NEMO/LIM_SRC/iceini.F90

    r712 r717  
    6363                  
    6464      ! Louvain la Neuve Ice model 
    65       IF( nacc == 1 ) THEN 
    66           dtsd2   = nfice * rdtmin * 0.5 
    67           rdt_ice = nfice * rdtmin 
    68       ELSE 
    69           dtsd2   = nfice * rdt * 0.5 
    70           rdt_ice = nfice * rdt 
    71       ENDIF 
     65      dtsd2   = nn_fsbc * rdttra(1) * 0.5 
     66      rdt_ice = nn_fsbc * rdttra(1) 
    7267 
    7368      CALL lim_msh                    ! ice mesh initialization 

  • trunk/NEMO/LIM_SRC/limdia.F90

    r699 r717  
    1919   USE par_ice         ! ice parameters 
    2020   USE ice_oce         ! ice variables 
     21   USE sbc_oce         ! surface boundary condition variables 
    2122   USE daymod          ! 
    2223   USE dom_ice         ! 
    … …  
    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/limdyn.F90

    r716 r717  
    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_ice_lim 
    712   !!---------------------------------------------------------------------- 
    … …  
    1116   !!    lim_dyn_init : initialization and namelist read 
    1217   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    14    USE phycst 
    15    USE in_out_manager  ! I/O manager 
    16    USE dom_ice 
    17    USE dom_oce         ! ocean space and time domain 
    18    USE ice 
    19    USE ice_oce 
    20    USE iceini 
    21    USE limistate 
    22    USE limrhg          ! 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            ! 
     22   USE ice_oce        ! 
     23   USE dom_ice        ! 
     24   USE iceini         ! 
     25   USE limistate      ! 
     26   USE limrhg         ! 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  ! routine called by ice_step 
    32  
    33    !! * Module variables 
     36   PUBLIC   lim_dyn   ! routine called by ice_step 
     37 
    3438   REAL(wp)  ::  rone    = 1.e0   ! constant value 
    3539 
    36  
    3740#  include "vectopt_loop_substitute.h90" 
    38  
    39    !!---------------------------------------------------------------------- 
    40    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
     41   !!---------------------------------------------------------------------- 
     42   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
    4143   !! $Id$ 
    42    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     44   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4345   !!---------------------------------------------------------------------- 
    4446 
    … …  
    5759      !!              - computation of the stress at the ocean surface          
    5860      !!              - treatment of the case if no ice dynamic 
    59       !! History : 
    60       !!   1.0  !  01-04  (LIM)  Original code 
    61       !!   2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp 
    6261      !!--------------------------------------------------------------------- 
    6362      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    64  
    65       INTEGER ::   ji, jj             ! dummy loop indices 
    66       INTEGER ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology 
    67       REAL(wp) ::   & 
    68          ztairx, ztairy,           &  ! tempory scalars 
    69          zsang , zrhomod,             & 
    70          ztglx , ztgly ,           & 
    71          zt11, zt12, zt21, zt22 ,  & 
    72          zustm, zsfrld, zsfrldm4,  & 
    73          zu_ice, zv_ice, ztair2 
    74       REAL(wp),DIMENSION(jpi,jpj) ::   zmod 
    75       REAL(wp),DIMENSION(jpj) ::   & 
    76          zind,                     &  ! i-averaged indicator of sea-ice 
    77          zmsk                         ! i-averaged of tmask 
     63      !! 
     64      INTEGER  ::   ji, jj             ! dummy loop indices 
     65      INTEGER  ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology 
     66      REAL(wp) ::   zcoef              ! temporary scalar 
     67      REAL(wp), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice 
     68      REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask 
     69      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io   ! ice-ocean velocity 
     70!!$      INTEGER ::   ji, jj             ! dummy loop indices 
     71!!$      INTEGER ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology 
     72!!$      REAL(wp) ::   & 
     73!!$         ztairx, ztairy,           &  ! tempory scalars 
     74!!$         zsang , zrhomod,             & 
     75!!$         ztglx , ztgly ,           & 
     76!!$         zt11, zt12, zt21, zt22 ,  & 
     77!!$         zustm, zsfrld, zsfrldm4,  & 
     78!!$         zu_ice, zv_ice, ztair2 
     79!!$      REAL(wp),DIMENSION(jpi,jpj) ::   zmod 
     80!!$      REAL(wp),DIMENSION(jpj) ::   & 
     81!!$         zind,                     &  ! i-averaged indicator of sea-ice 
     82!!$         zmsk                         ! i-averaged of tmask 
    7883      !!--------------------------------------------------------------------- 
    7984 
    80       IF( kt == nit000  )   CALL lim_dyn_init   ! Initialization (first time-step only) 
     85      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only) 
    8186       
    82       IF ( ln_limdyn ) THEN 
     87      IF( ln_limdyn ) THEN 
    8388 
    8489         ! Mean ice and snow thicknesses.           
    … …  
    8691         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:) 
    8792 
    88          u_io(:,:)  = u_io(:,:) * tmu(:,:) 
    89          v_io(:,:)  = v_io(:,:) * tmu(:,:) 
    90         
    9193         !                                         ! Rheology (ice dynamics) 
    9294         !                                         ! ======== 
    … …  
    98100            i_j1 = 1    
    99101            i_jpj = jpj 
    100             IF(ln_ctl)    THEN 
    101                CALL prt_ctl_info('lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj) 
    102             ENDIF 
     102            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    103103            CALL lim_rhg( i_j1, i_jpj ) 
    104104 
    … …  
    108108               zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line 
    109109               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line 
    110    !!i         write(numout,*) narea, 'limdyn' , jj, zind(jj), zmsk(jj) 
    111110            END DO 
    112111 
    … …  
    159158         ENDIF 
    160159 
    161          IF(ln_ctl)   THEN  
    162             CALL prt_ctl(tab2d_1=u_io , clinfo1=' lim_dyn  : u_io :', tab2d_2=v_io , clinfo2=' v_io :') 
    163             CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_dyn  : u_ice:', tab2d_2=v_ice, clinfo2=' v_ice:') 
    164          ENDIF 
    165           
    166          !                                         ! Ice-Ocean stress 
    167          !                                         ! ================         
    168          DO jj = 1, jpj 
    169             DO ji = 1, jpi 
    170 !!              zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg    ! do the full loop and avoid lbc_lnk 
    171                zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg 
    172                zu_ice   = u_ice(ji,jj) - u_io(ji,jj) 
    173                zv_ice   = v_ice(ji,jj) - v_io(ji,jj) 
    174                zrhomod  = zu_ice * zu_ice + zv_ice * zv_ice  
    175                zmod (ji,jj) = zrhomod  
    176                zrhomod  = rhoco * SQRT( zrhomod )  
    177                ftaux(ji,jj) = zrhomod * ( cangvg * zu_ice - zsang * zv_ice )  
    178                ftauy(ji,jj) = zrhomod * ( cangvg * zv_ice + zsang * zu_ice )  
     160         IF(ln_ctl)   CALL prt_ctl(tab2d_1=ui_ice , clinfo1=' lim_dyn  : ui_ice :', tab2d_2=vi_ice , clinfo2=' vi_ice :') 
     161          
     162!!$         !                                         ! Ice-Ocean stress 
     163!!$         !                                         ! ================         
     164!!$         DO jj = 1, jpj 
     165!!$            DO ji = 1, jpi 
     166!!$!!              zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg    ! do the full loop and avoid lbc_lnk 
     167!!$               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg 
     168!!$               zu_ice   = u_ice(ji,jj) - u_io(ji,jj) 
     169!!$               zv_ice   = v_ice(ji,jj) - v_io(ji,jj) 
     170!!$               zrhomod  = zu_ice * zu_ice + zv_ice * zv_ice  
     171!!$               zmod (ji,jj) = zrhomod  
     172!!$               zrhomod  = rhoco * SQRT( zrhomod )  
     173!!$               ftaux(ji,jj) = zrhomod * ( cangvg * zu_ice - zsang * zv_ice )  
     174!!$               ftauy(ji,jj) = zrhomod * ( cangvg * zv_ice + zsang * zu_ice )  
     175!!$            END DO 
     176!!$         END DO 
     177 
     178         ! computation of friction velocity 
     179         ! -------------------------------- 
     180         ! ice-ocean velocity at U & V-points (ui_ice vi_ice at I-point ; ssu_m, ssv_m at U- & V-points) 
     181          
     182         DO jj = 1, jpjm1 
     183            DO ji = 1, fs_jpim1   ! vector opt. 
     184               zu_io(ji,jj) = 0.5 * ( ui_ice(ji+1,jj+1) + ui_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
     185               zv_io(ji,jj) = 0.5 * ( vi_ice(ji+1,jj+1) + vi_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
    179186            END DO 
    180187         END DO 
    181  
    182          ! computation of friction velocity 
     188         ! frictional velocity at T-point 
    183189         DO jj = 2, jpjm1 
    184             DO ji = fs_2, fs_jpim1     ! vector opt. 
    185                ust2s(ji,jj) = 0.25 * cw * ( zmod(ji,jj+1) + zmod(ji+1,jj+1) +    & 
    186                   &                         zmod(ji,jj  ) + zmod(ji+1,jj  ) ) * tms(ji,jj) 
     190            DO ji = fs_2, fs_jpim1   ! vector opt. 
     191               ust2s(ji,jj) = 0.5 * cw                                                          & 
     192                  &         * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
     193                  &            + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj) 
    187194            END DO 
    188195         END DO 
    189  
    190        ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
    191                      
    192           ftaux(:,:) = gtaux(:,:)  
    193           ftauy(:,:) = gtauy(:,:)  
    194  
    195           DO jj = 2, jpjm1 
    196              DO ji = 2, jpim1 
    197                 ztair2 = gtaux(ji  ,jj+1) * gtaux(ji  ,jj+1) + gtauy(ji  ,jj+1) * gtauy(ji  ,jj+1)   & 
    198                 &      + gtaux(ji+1,jj+1) * gtaux(ji+1,jj+1) + gtauy(ji+1,jj+1) * gtauy(ji+1,jj+1)   & 
    199                 &      + gtaux(ji  ,jj  ) * gtaux(ji  ,jj  ) + gtauy(ji  ,jj  ) * gtauy(ji  ,jj  )   & 
    200                 &      + gtaux(ji+1,jj  ) * gtaux(ji+1,jj  ) + gtauy(ji+1,jj  ) * gtauy(ji+1,jj  )  
    201  
    202                 ust2s(ji,jj) = 0.25 / rau0 * SQRT( ztair2 ) * tms(ji,jj) 
     196!!$         DO jj = 2, jpjm1 
     197!!$            DO ji = fs_2, fs_jpim1     ! vector opt. 
     198!!$               ust2s(ji,jj) = 0.25 * cw * ( zmod(ji,jj+1) + zmod(ji+1,jj+1) +    & 
     199!!$                  &                         zmod(ji,jj  ) + zmod(ji+1,jj  ) ) * tms(ji,jj) 
     200!!$            END DO 
     201!!$         END DO 
     202         ! 
     203      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
     204         ! 
     205         zcoef = SQRT( 0.5 ) / rau0 
     206         DO jj = 2, jpjm1 
     207            DO ji = fs_2, fs_jpim1   ! vector opt. 
     208               ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
     209                  &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) 
    203210            END DO 
    204211         END DO 
    205  
     212         ! 
    206213      ENDIF 
    207214 
    208215      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point 
    209216 
    210       IF(ln_ctl) THEN  
    211             CALL prt_ctl(tab2d_1=ftaux , clinfo1=' lim_dyn  : ftaux :', tab2d_2=ftauy , clinfo2=' ftauy :') 
    212             CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    213       ENDIF 
     217      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    214218 
    215219   END SUBROUTINE lim_dyn 
    … …  
    220224      !!                  ***  ROUTINE lim_dyn_init  *** 
    221225      !! 
    222       !! ** Purpose : Physical constants and parameters linked to the ice 
    223       !!      dynamics 
    224       !! 
    225       !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic 
    226       !!       parameter values called at the first timestep (nit000) 
     226      !! ** Purpose :   Physical constants and parameters linked to the ice 
     227      !!              dynamics 
     228      !! 
     229      !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic 
     230      !!              parameter values 
    227231      !! 
    228232      !! ** input   :   Namelist namicedyn 
    229       !! 
    230       !! history : 
    231       !!  8.5  ! 03-08 (C. Ethe) original code 
    232233      !!------------------------------------------------------------------- 
    233234      NAMELIST/namicedyn/ epsd, alpha,     & 
    … …  
    236237      !!------------------------------------------------------------------- 
    237238 
    238       ! Define the initial parameters 
    239       ! ------------------------- 
    240  
    241       ! Read Namelist namicedyn 
    242       REWIND ( numnam_ice ) 
     239      REWIND ( numnam_ice )                       ! Read Namelist namicedyn 
    243240      READ   ( numnam_ice  , namicedyn ) 
    244       IF(lwp) THEN 
     241 
     242      IF(lwp) THEN                                ! Control print 
    245243         WRITE(numout,*) 
    246244         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 
    … …  
    254252         WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl 
    255253         WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw 
    256          WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg 
     254         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg, ' degrees' 
    257255         WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar 
    258256         WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg 
    … …  
    263261      ENDIF 
    264262 
    265       !  Initialization 
    266263      usecc2 = 1.0 / ( ecc * ecc ) 
    267264      rhoco  = rau0 * cw 
    268       angvg  = angvg * rad 
     265      angvg  = angvg * rad      ! convert angvg from degree to radian 
    269266      sangvg = SIN( angvg ) 
    270267      cangvg = COS( angvg ) 
    271268      pstarh = pstar / 2.0 
    272       sdvt(:,:) = 0.e0 
    273  
    274       !  Diffusion coefficients. 
    275       ahiu(:,:) = ahi0 * umask(:,:,1) 
     269      ! 
     270      ahiu(:,:) = ahi0 * umask(:,:,1)            ! Ice eddy Diffusivity coefficients. 
    276271      ahiv(:,:) = ahi0 * vmask(:,:,1) 
    277  
     272      ! 
    278273   END SUBROUTINE lim_dyn_init 
    279274 

  • trunk/NEMO/LIM_SRC/limistate.F90

    r699 r717  
    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_ice_lim 
    … …  
    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         ! ice parameters 
    2323   USE ice_oce         ! ice variables 
    2424   USE dom_ice 
    25    USE ice             ! ??? 
    2625   USE lbclnk 
     26   USE oce 
    2727   USE ice 
    2828   USE iom 
    … …  
    3232   PRIVATE 
    3333 
    34    PUBLIC lim_istate      ! routine called by lim_init.F90 
    35  
    36    REAL(wp) ::           &  !!! ** init namelist (namiceini) ** 
     34   PUBLIC   lim_istate   ! routine called by lim_init.F90 
     35 
     36   REAL(wp) ::           & !!! ** init namelist (namiceini) ** 
    3737      ttest  = 2.0  ,    &  ! threshold water temperature for initial sea ice 
    3838      hninn  = 0.5  ,    &  ! initial snow thickness in the north 
    … …  
    6767      REAL(wp), DIMENSION(jpi,jpj) ::   ztn   ! workspace 
    6868      !-------------------------------------------------------------------- 
    69  
    70        CALL lim_istate_init     !  reading the initials parameters of the ice 
    71  
    72       !-- Initialisation of sst,sss,u,v do i=1,jpi 
    73       u_io(:,:)  = 0.e0       ! ice velocity in x direction 
    74       v_io(:,:)  = 0.e0       ! ice velocity in y direction 
    75  
    76       IF( ln_limini ) THEN    !  
    77          
    78          ! Initialisation at tn if no ice or sst_ini if ice 
    79          ! Idem for salinity 
    80  
    81       !--- Criterion for presence (zidto=1.) or absence (zidto=0.) of ice 
    82          DO jj = 1 , jpj 
    83             DO ji = 1 , jpi 
    84                 
    85                zidto = MAX(zzero, - SIGN(1.,frld(ji,jj) - 1.)) 
    86                 
    87                sst_io(ji,jj) = ( nfice - 1 ) * (zidto * sst_ini(ji,jj)  + &   ! use the ocean initial values 
    88                     &          (1.0 - zidto ) * ( tn(ji,jj,1) + rt0 ))        ! tricky trick *(nfice-1) ! 
    89                sss_io(ji,jj) = ( nfice - 1 ) * (zidto * sss_ini(ji,jj) + & 
    90                     &          (1.0 - zidto ) *  sn(ji,jj,1) ) 
    91  
    92                ! to avoid the the melting of ice, several layers (mixed layer) should be 
    93                ! set to sst_ini (sss_ini) if there is ice 
    94                ! example for one layer  
    95                ! tn(ji,jj,1) = zidto * ( sst_ini(ji,jj) - rt0 )  + (1.0 - zidto ) *  tn(ji,jj,1) 
    96                ! sn(ji,jj,1) = zidto * sss_ini(ji,jj)  + (1.0 - zidto ) *  sn(ji,jj,1) 
    97                ! tb(ji,jj,1) = tn(ji,jj,1) 
    98                ! sb(ji,jj,1) = sn(ji,jj,1) 
    99             END DO 
    100          END DO 
    101           
    102           
    103          !  tfu: Melting point of sea water 
    104          tfu(:,:)  = ztf    
    105           
    106          tfu(:,:)  = ABS ( rt0 - 0.0575       * sss_ini(:,:)                               & 
    107               &                    + 1.710523e-03 * sss_ini(:,:) * SQRT( sss_ini(:,:) )    & 
    108               &                    - 2.154996e-04 * sss_ini(:,:) * sss_ini(:,:) ) 
    109       ELSE                     ! 
    110  
     69  
     70      CALL lim_istate_init     !  reading the initials parameters of the ice 
     71 
     72      IF( .NOT. ln_limini ) THEN   
    11173          
    11274         ! Initialisation at tn or -2 if ice 
    … …  
    11779            END DO 
    11880         END DO 
    119           
    120          u_io  (:,:) = 0.e0 
    121          v_io  (:,:) = 0.e0 
    122          sst_io(:,:) = ( nfice - 1 ) * ( tn(:,:,1) + rt0 )   ! use the ocean initial values 
    123          sss_io(:,:) = ( nfice - 1 ) *   sn(:,:,1)           ! tricky trick *(nfice-1) ! 
    124           
    125          ! reference salinity 34psu 
     81                   
     82         !  tfu: Melting point of sea water [Kelvin] 
    12683         zs0 = 34.e0 
    127          ztf = ABS ( rt0 - 0.0575       * zs0                           & 
    128               &                    + 1.710523e-03 * zs0 * SQRT( zs0 )   & 
    129               &                    - 2.154996e-04 * zs0 *zs0          ) 
    130           
    131          !  tfu: Melting point of sea water 
    132          tfu(:,:)  = ztf    
     84         ztf = rt0 + ( - 0.0575 + 1.710523e-3 * SQRT( zs0 ) - 2.154996e-4 * zs0 ) * zs0 
     85         tfu(:,:) = ztf 
    13386          
    13487         DO jj = 1, jpj 
    … …  
    153106         tbif  (:,:,2) = tfu(:,:) 
    154107         tbif  (:,:,3) = tfu(:,:) 
    155        
     108 
    156109      ENDIF 
     110      
    157111      fsbbq (:,:)   = 0.e0 
    158112      qstoif(:,:)   = 0.e0 
    159       u_ice (:,:)   = 0.e0 
    160       v_ice (:,:)   = 0.e0 
     113      ui_ice(:,:)   = 0.e0 
     114      vi_ice(:,:)   = 0.e0 
    161115# if defined key_coupled 
    162116      albege(:,:)   = 0.8 * tms(:,:) 
    … …  
    192146 
    193147      CALL lbc_lnk( hsnif, 'T', 1. ) 
    194       CALL lbc_lnk( sist , 'T', 1. ) 
     148      CALL lbc_lnk( sist , 'T', 1. , pval = rt0 )      ! set rt0 on closed boundary (required by bulk formulation) 
    195149      DO jk = 1, jplayersp1 
    196150         CALL lbc_lnk(tbif(:,:,jk), 'T', 1. ) 
    … …  
    198152      CALL lbc_lnk( fsbbq  , 'T', 1. ) 
    199153      CALL lbc_lnk( qstoif , 'T', 1. ) 
    200       CALL lbc_lnk( sss_io , 'T', 1. ) 
    201       ! 
     154 
    202155   END SUBROUTINE lim_istate 
    203156 
    … …  
    210163      !! 
    211164      !! ** Method  :   Read the namiceini namelist and check the parameter  
    212       !!                values called at the first timestep (nit000) 
    213       !!                or 
    214       !!                Read 7 variables from a previous restart file 
    215       !!                sst, sst, hicif, hsnif, frld, ts & tbif 
     165      !!       values called at the first timestep (nit000) 
    216166      !! 
    217167      !! ** input   :   Namelist namiceini 
    … …  
    223173         &                hnins, hgins, alins 
    224174      !!------------------------------------------------------------------- 
    225        
    226       ! Read Namelist namiceini  
    227       REWIND ( numnam_ice ) 
     175      ! 
     176      REWIND ( numnam_ice )               ! Read Namelist namiceini  
    228177      READ   ( numnam_ice , namiceini ) 
    229        
    230       IF(.NOT. ln_limini) THEN  
    231          IF(lwp) THEN 
    232             WRITE(numout,*) 
    233             WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 
    234             WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    235             WRITE(numout,*) '         threshold water temp. for initial sea-ice    ttest      = ', ttest 
    236             WRITE(numout,*) '         initial snow thickness in the north          hninn      = ', hninn 
    237             WRITE(numout,*) '         initial ice thickness in the north           hginn      = ', hginn  
    238             WRITE(numout,*) '         initial leads area in the north              alinn      = ', alinn             
    239             WRITE(numout,*) '         initial snow thickness in the south          hnins      = ', hnins  
    240             WRITE(numout,*) '         initial ice thickness in the south           hgins      = ', hgins 
    241             WRITE(numout,*) '         initial leads area in the south              alins      = ', alins 
    242          ENDIF 
     178 
     179      IF(lwp) THEN 
     180         WRITE(numout,*) 
     181         WRITE(numout,*) 'lim_istate_init : 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 
    243191      ENDIF 
    244192 
    245193      IF( ln_limini ) THEN                      ! Ice initialization using input file 
    246  
     194         ! 
    247195         CALL iom_open( 'Ice_initialization.nc', inum_ice ) 
    248  
     196         ! 
    249197         IF( inum_ice > 0 ) THEN 
    250             IF(lwp) THEN 
    251                WRITE(numout,*) ' ' 
    252                WRITE(numout,*) 'lim_istate_init : ice state initialization with : Ice_initialization.nc' 
    253                WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    254                WRITE(numout,*) '         Ice state initialization using input file    ln_limini  = ', ln_limini 
    255                WRITE(numout,*) ' ' 
    256             ENDIF 
     198            IF(lwp) WRITE(numout,*) 
     199            IF(lwp) WRITE(numout,*) '                  ice state initialization with : Ice_initialization.nc' 
    257200             
    258             CALL iom_get( inum_ice, jpdom_data, 'sst'  , sst_ini(:,:) )         
    259             CALL iom_get( inum_ice, jpdom_data, 'sss'  , sss_ini(:,:) )        
    260             CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif  (:,:) )       
    261             CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif  (:,:) )       
    262             CALL iom_get( inum_ice, jpdom_data, 'frld' , frld   (:,:) )      
    263             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  ) 
    264205            CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif(1:nlci,1:nlcj,:),   & 
    265206                 &        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) ) 
    … …  
    269210 
    270211            CALL iom_close( inum_ice) 
    271              
     212            ! 
    272213         ENDIF 
    273214      ENDIF 
    274       ! 
     215      !      
    275216   END SUBROUTINE lim_istate_init 
    276217 

  • trunk/NEMO/LIM_SRC/limrhg.F90

    r701 r717  
    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    !!             3.0  !  07-06  (S. Masson, G. Madec)  ice/atmosphere stress 
    10    !!                             provided at I-point in forced and coupled mode 
     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 
    1111   !!---------------------------------------------------------------------- 
    1212#if defined key_ice_lim 
    … …  
    1414   !!   'key_ice_lim'                                     LIM sea-ice model 
    1515   !!---------------------------------------------------------------------- 
     16   !!---------------------------------------------------------------------- 
    1617   !!   lim_rhg   : computes ice velocities 
    1718   !!---------------------------------------------------------------------- 
    18    !! * Modules used 
    19    USE phycst 
    20    USE par_oce 
    21    USE ice_oce         ! ice variables 
    22    USE dom_ice 
    23    USE ice 
    24    USE lbclnk 
    25    USE lib_mpp 
    26    USE in_out_manager  ! I/O manager 
    27    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        ! domaine: ice variables 
     23   USE phycst         ! physical constant 
     24   USE ice            ! 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 
    2829 
    2930   IMPLICIT NONE 
    3031   PRIVATE 
    3132 
    32    !! * Routine accessibility 
    33    PUBLIC lim_rhg  ! routine called by lim_dyn 
    34  
    35    !! * Module variables 
    36    REAL(wp)  ::           &  ! constant values 
    37       rzero   = 0.e0   ,  & 
    38       rone    = 1.e0 
    39    !!---------------------------------------------------------------------- 
    40    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    41    !! $Id$  
    42    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     33   PUBLIC   lim_rhg   ! 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   !! $Header: $ 
     43   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4344   !!---------------------------------------------------------------------- 
    4445 
    … …  
    5455      !!  viscous-plastic law including shear strength and a bulk rheology. 
    5556      !! 
    56       !! ** 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 
    5762      !! 
     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 
    5891      !!------------------------------------------------------------------- 
    59       INTEGER, INTENT(in) ::   k_j1 ,  &  ! southern j-index for ice computation 
    60          &                     k_jpj      ! northern j-index for ice computation 
    61  
    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, zic 
    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(:,:) 
    93  
    94       ! masked Ice mass, ice strength, Ice-cover fraction and Lead faction  at T-point  
     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 
     132 
     133      ! Ice mass, ice strength, and wind stress at the center            | 
     134      ! of the grid squares.                                             | 
    95135      !------------------------------------------------------------------- 
    96136 
     137!CDIR NOVERRCHK 
    97138      DO jj = k_j1 , k_jpj-1 
     139!CDIR NOVERRCHK 
    98140         DO ji = 1 , jpi 
    99             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) ) 
    100145            zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) ) 
    101  
    102             zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )     
    103             zb2(ji,jj)    = tms(ji,jj) *         frld(ji,jj)       
     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) ) 
     149 
     150!!gm   #if defined key_lim_cp1 && defined key_coupled 
     151!!gm        zi1(ji,jj)    = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) ) 
     152!!gm        zi2(ji,jj)    = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) ) 
     153!!gm   #else 
     154!!gm        zi1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 
     155!!gm        zi2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 
     156!!gm   #endif 
    104157         END DO 
    105158      END DO 
    106159 
    107       !  Wind stress, coriolis, mass and Gradient of ice strenght at I-point 
     160 
     161      !--------------------------------------------------------------------------- 
     162      !  Wind stress, coriolis and mass terms at the corners of the grid squares | 
     163      !  Gradient of ice strenght.                                               | 
    108164      !--------------------------------------------------------------------------- 
    109165          
    110166      DO jj = k_j1+1, k_jpj-1 
    111          DO ji = 2, jpi 
    112             zstms = tms(ji,jj  ) * wght(ji,jj,2,2) + tms(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    113                &  + tms(ji,jj-1) * wght(ji,jj,2,1) + tms(ji-1,jj-1) * wght(ji,jj,1,1) 
     167         DO ji = fs_2, jpi 
     168            zstms = zztms(ji,jj  ) * wght(ji,jj,2,2) + zztms(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     169               &  + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1) 
    114170            zusw  = 1.0 / MAX( zstms, epsd ) 
    115171 
    116             ! Leads area at I-point 
    117             zfrld(ji,jj) = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    118                &           + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
    119  
    120             ! Ice cover area at I-point 
    121             zic          = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    122                &           + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
     172            zt11 = zztms(ji  ,jj  ) * zzfrld(ji  ,jj  )  
     173            zt12 = zztms(ji-1,jj  ) * zzfrld(ji-1,jj  )  
     174            zt21 = zztms(ji  ,jj-1) * zzfrld(ji  ,jj-1)  
     175            zt22 = zztms(ji-1,jj-1) * zzfrld(ji-1,jj-1) 
     176 
     177            ! Leads area. 
     178            zfrld(ji,jj) =  (  zt11 * wght(ji,jj,2,2) + zt12 * wght(ji,jj,1,2)   & 
     179               &             + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw 
     180 
     181            ! Mass and coriolis coeff. at I-point 
     182            zmass(ji,jj) = ( zmasst(ji,jj  ) * wght(ji,jj,2,2) + zmasst(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     183               &           + zmasst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
     184            zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj) 
    123185 
    124186            ! Wind stress. 
    125             ztagnx = zic * gtaux(ji,jj) 
    126             ztagny = zic * gtauy(ji,jj) 
    127  
    128             ! Mass and coriolis coeff. 
    129             zmass(ji,jj) = ( za1(ji,jj  ) * wght(ji,jj,2,2) + za1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    130                &           + za1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
    131             zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj) 
     187!!gm   #if defined key_lim_cp1 && defined key_coupled 
     188!!gm   ztagnx = ( zi1(ji,jj  ) * wght(ji,jj,2,2) + zi1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     189!!gm      &     + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
     190!!gm   ztagny = ( zi2(ji,jj  ) * wght(ji,jj,2,2) + zi2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     191!!gm      &     + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 
     192!!gm   #else 
     193!!gm   ztagnx = ( zi1(ji,jj  ) * wght(ji,jj,2,2) + zi1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     194!!gm      &     + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj) 
     195!!gm   ztagny = ( zi2(ji,jj  ) * wght(ji,jj,2,2) + zi2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     196!!gm      &     + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj) 
     197!!gm   #endif 
     198            !!gm  always stress provided at I-point (ocean F-point) 
     199            ztagnx = ( zi1(ji,jj  ) * wght(ji,jj,2,2) + zi1(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     200               &     + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utaui_ice(ji,jj) 
     201            ztagny = ( zi2(ji,jj  ) * wght(ji,jj,2,2) + zi2(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     202               &     + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtaui_ice(ji,jj) 
    132203 
    133204            ! Gradient of ice strength 
    … …  
    143214 
    144215            ! Computation of the velocity field taking into account the ice-ice interaction.                                  
    145             ! Terms that are independent of the velocity field. 
    146             za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_io(ji,jj) - zgphsx 
    147             za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_io(ji,jj) - zgphsy 
     216            ! Terms that are independent of the ice velocity field. 
     217            za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * vi_oce(ji,jj) - zgphsx 
     218            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * ui_oce(ji,jj) - zgphsy 
    148219         END DO 
    149220      END DO 
    150  
    151 !! inutile!! 
    152 !!??    CALL lbc_lnk( za1ct, 'I', -1. ) 
    153 !!??    CALL lbc_lnk( za2ct, 'I', -1. ) 
    154221 
    155222 
    … …  
    164231         ! Computation of free drift field for free slip boundary conditions. 
    165232 
    166            DO jj = k_j1, k_jpj-1 
    167               DO ji = 1, jpim1 
    168                  !- Rate of strain tensor. 
    169                  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) )  & 
    170                     &   + 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) ) 
    171                  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) )  & 
    172                     &   - 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) ) 
    173                  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) )  & 
    174                     &   + 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) ) 
    175                  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) )  & 
    176                     &   - 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) ) 
    177  
    178                  !- Rate of strain tensor.  
    179                  zdgp = zt11 + zt22 
    180                  zdgi = zt12 + zt21 
    181                  ztrace2 = zdgp * zdgp  
    182                  zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi 
    183  
    184                  !  Creep limit depends on the size of the grid. 
    185                  zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2),  creepl) 
    186  
    187                  !-  Computation of viscosities. 
    188                  zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 
    189                  zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 
    190               END DO 
    191            END DO 
    192 !!??       CALL lbc_lnk( zviszeta, 'I', -1. )  ! or T point???   semble reellement inutile 
    193 !!??       CALL lbc_lnk( zviseta , 'I', -1. ) 
    194  
    195  
    196            !-  Determination of zc1u, zc2u, zc1v and zc2v. 
    197            DO jj = k_j1+1, k_jpj-1 
    198               DO ji = 2, jpim1 
    199                  ze11   =  akappa(ji-1,jj-1,1,1) 
    200                  ze12   = +akappa(ji-1,jj-1,2,2) 
    201                  ze22   =  akappa(ji-1,jj-1,2,1) 
    202                  ze21   = -akappa(ji-1,jj-1,1,2) 
    203                  zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
    204                  zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
    205                  zvis12 =       zviseta (ji-1,jj-1) + dm 
    206                  zvis21 =       zviseta (ji-1,jj-1) 
    207  
    208                  zdiag = zvis22 * ( ze11 + ze22 ) 
    209                  zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag 
    210                  zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21 
    211                  zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag 
    212                  zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12 
    213  
    214                  ze11   = -akappa(ji,jj-1,1,1) 
    215                  ze12   = +akappa(ji,jj-1,2,2) 
    216                  ze22   =  akappa(ji,jj-1,2,1) 
    217                  ze21   = -akappa(ji,jj-1,1,2) 
    218                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
    219                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
    220                  zvis12 =       zviseta (ji,jj-1) + dm 
    221                  zvis21 =       zviseta (ji,jj-1) 
    222  
    223                  zdiag = zvis22 * ( ze11 + ze22 ) 
    224                  zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag 
    225                  zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21 
    226                  zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag 
    227                  zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12 
    228  
    229                  ze11   =  akappa(ji-1,jj,1,1) 
    230                  ze12   = -akappa(ji-1,jj,2,2) 
    231                  ze22   =  akappa(ji-1,jj,2,1) 
    232                  ze21   = -akappa(ji-1,jj,1,2) 
    233                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
    234                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
    235                  zvis12 =       zviseta (ji-1,jj) + dm 
    236                  zvis21 =       zviseta (ji-1,jj) 
    237  
    238                  zdiag = zvis22 * ( ze11 + ze22 )  
    239                  zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag 
    240                  zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21 
    241                  zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag 
    242                  zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12 
    243  
    244                  ze11   = -akappa(ji,jj,1,1) 
    245                  ze12   = -akappa(ji,jj,2,2) 
    246                  ze22   =  akappa(ji,jj,2,1) 
    247                  ze21   = -akappa(ji,jj,1,2) 
    248                  zvis11 = 2.0 * zviseta (ji,jj) + dm 
    249                  zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
    250                  zvis12 =       zviseta (ji,jj) + dm 
    251                  zvis21 =       zviseta (ji,jj) 
    252  
    253                  zdiag = zvis22 * ( ze11 + ze22 ) 
    254                  zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag 
    255                  zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21 
    256                  zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag  
    257                  zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12 
    258               END DO 
    259            END DO 
    260  
    261            DO jj = k_j1+1, k_jpj-1 
    262               DO ji = 2, jpim1 
    263                  zc1u(ji,jj) =   & 
    264                     + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   & 
    265                     - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   & 
    266                     - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   & 
    267                     + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   & 
    268                     + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   & 
    269                     + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   & 
    270                     - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   & 
    271                     - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 
    272                   
    273                  zc2u(ji,jj) =   & 
    274                     + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   & 
    275                     - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   & 
    276                     - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   & 
    277                     + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   & 
    278                     - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   & 
    279                     - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   & 
    280                     + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   & 
    281                     + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 
    282              END DO 
    283            END DO 
    284  
    285            DO jj = k_j1+1, k_jpj-1 
    286               DO ji = 2, jpim1 
    287                  !  zc1v , zc2v. 
    288                  ze11   =  akappa(ji-1,jj-1,1,2) 
    289                  ze12   = -akappa(ji-1,jj-1,2,1) 
    290                  ze22   = +akappa(ji-1,jj-1,2,2) 
    291                  ze21   =  akappa(ji-1,jj-1,1,1) 
    292                  zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
    293                  zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
    294                  zvis12 =       zviseta (ji-1,jj-1) + dm 
    295                  zvis21 =       zviseta (ji-1,jj-1) 
    296  
    297                  zdiag = zvis22 * ( ze11 + ze22 ) 
    298                  zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag 
    299                  zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21 
    300                  zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag 
    301                  zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12 
     233!CDIR NOVERRCHK 
     234         DO jj = k_j1, k_jpj-1 
     235!CDIR NOVERRCHK 
     236            DO ji = 1, fs_jpim1 
     237               !- Rate of strain tensor. 
     238               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) )  & 
     239                  &   + 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) ) 
     240               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) )  & 
     241                  &   - 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) ) 
     242               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) )  & 
     243                  &   + 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) ) 
     244               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) )  & 
     245                  &   - 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) ) 
     246 
     247               !- Rate of strain tensor.  
     248               zdgp = zt11 + zt22 
     249               zdgi = zt12 + zt21 
     250               ztrace2 = zdgp * zdgp  
     251               zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi 
     252 
     253               !  Creep limit depends on the size of the grid. 
     254               zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2 ),  creepl) 
     255 
     256               !-  Computation of viscosities. 
     257               zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 
     258               zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 
     259            END DO 
     260         END DO 
     261 
     262         !-  Determination of zc1u, zc2u, zc1v and zc2v. 
     263         DO jj = k_j1+1, k_jpj-1 
     264            DO ji = fs_2, fs_jpim1 
     265               !* zc1u , zc2v 
     266               zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
     267               zvis12 =       zviseta (ji-1,jj-1) + dm 
     268               zvis21 =       zviseta (ji-1,jj-1) 
     269               zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
     270               zdiag  = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) ) 
     271               zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,1) + zdiag 
     272               zs12_11 =  zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2) 
     273               zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2) 
     274               zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,1) + zdiag 
     275 
     276               zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
     277               zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
     278               zvis12 =       zviseta (ji,jj-1) + dm 
     279               zvis21 =       zviseta (ji,jj-1) 
     280               zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) ) 
     281               zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag 
     282               zs12_21 =  zvis12 * akappa(ji,jj-1,2,2) - zvis21 * akappa(ji,jj-1,1,2) 
     283               zs22_21 =  zvis11 * akappa(ji,jj-1,2,1) + zdiag 
     284               zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2) 
     285 
     286               zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
     287               zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
     288               zvis12 =       zviseta (ji-1,jj) + dm 
     289               zvis21 =       zviseta (ji-1,jj) 
     290               zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) ) 
     291               zs11_12 =  zvis11 * akappa(ji-1,jj,1,1) + zdiag 
     292               zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2) 
     293               zs22_12 =  zvis11 * akappa(ji-1,jj,2,1) + zdiag 
     294               zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2) 
     295 
     296               zvis11 = 2.0 * zviseta (ji,jj) + dm 
     297               zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
     298               zvis12 =       zviseta (ji,jj) + dm 
     299               zvis21 =       zviseta (ji,jj) 
     300               zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) ) 
     301               zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag 
     302               zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2) 
     303               zs22_22 =  zvis11 * akappa(ji,jj,2,1) + zdiag 
     304               zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2) 
     305 
     306               zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22   & 
     307                  &          - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12   & 
     308                  &          - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11   & 
     309                  &          + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12   & 
     310                  &          + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21   & 
     311                  &          + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22   & 
     312                  &          - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21   & 
     313                  &          - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 
     314 
     315               zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22   & 
     316                  &          - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12   & 
     317                  &          - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11   & 
     318                  &          + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12   & 
     319                  &          - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21   & 
     320                  &          - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22   & 
     321                  &          + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21   & 
     322                  &          + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 
     323 
     324               !* zc1v , zc2v. 
     325               zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
     326               zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
     327               zvis12 =       zviseta (ji-1,jj-1) + dm 
     328               zvis21 =       zviseta (ji-1,jj-1) 
     329               zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) ) 
     330               zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,2) + zdiag 
     331               zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1) 
     332               zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,2) + zdiag 
     333               zs21_11 =  zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1) 
    302334  
    303                  ze11   =  akappa(ji,jj-1,1,2) 
    304                  ze12   = -akappa(ji,jj-1,2,1) 
    305                  ze22   = +akappa(ji,jj-1,2,2) 
    306                  ze21   = -akappa(ji,jj-1,1,1) 
    307                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
    308                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
    309                  zvis12 =       zviseta (ji,jj-1) + dm 
    310                  zvis21 =       zviseta (ji,jj-1) 
    311  
    312                  zdiag = zvis22 * ( ze11 + ze22 ) 
    313                  zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag 
    314                  zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21 
    315                  zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag 
    316                  zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12 
    317  
    318                  ze11   =  akappa(ji-1,jj,1,2) 
    319                  ze12   = -akappa(ji-1,jj,2,1) 
    320                  ze22   = -akappa(ji-1,jj,2,2) 
    321                  ze21   =  akappa(ji-1,jj,1,1) 
    322                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
    323                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
    324                  zvis12 =       zviseta (ji-1,jj) + dm 
    325                  zvis21 =       zviseta (ji-1,jj) 
    326  
    327                  zdiag = zvis22 * ( ze11 + ze22 ) 
    328                  zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag 
    329                  zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21 
    330                  zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag 
    331                  zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12 
    332  
    333                  ze11   =  akappa(ji,jj,1,2) 
    334                  ze12   = -akappa(ji,jj,2,1) 
    335                  ze22   = -akappa(ji,jj,2,2) 
    336                  ze21   = -akappa(ji,jj,1,1) 
    337                  zvis11 = 2.0 * zviseta (ji,jj) + dm 
    338                  zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
    339                  zvis12 =       zviseta (ji,jj) + dm 
    340                  zvis21 =       zviseta (ji,jj) 
    341  
    342                  zdiag = zvis22 * ( ze11 + ze22 ) 
    343                  zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag 
    344                  zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21 
    345                  zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag 
    346                  zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12 
    347  
    348               END DO 
    349            END DO 
    350  
    351            DO jj = k_j1+1, k_jpj-1 
    352               DO ji = 2, jpim1 
    353                  zc1v(ji,jj) =   & 
    354                     + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   & 
    355                     - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   & 
    356                     - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   & 
    357                     + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   & 
    358                     + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   & 
    359                     + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   & 
    360                     - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   & 
    361                     - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 
    362                  zc2v(ji,jj) =   & 
    363                     + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   & 
    364                     - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   & 
    365                     - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   & 
    366                     + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   & 
    367                     - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   & 
    368                     - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   & 
    369                     + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   & 
    370                     + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 
    371               END DO 
    372            END DO 
    373  
    374          ! Relaxation. 
    375             
    376 iflag:   DO jter = 1 , nbitdr 
    377  
    378             !  Store previous drift field.    
    379             DO jj = k_j1, k_jpj-1 
    380                zu_ice(:,jj) = u_ice(:,jj) 
    381                zv_ice(:,jj) = v_ice(:,jj) 
     335               zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
     336               zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
     337               zvis12 =       zviseta (ji,jj-1) + dm 
     338               zvis21 =       zviseta (ji,jj-1) 
     339               zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) ) 
     340               zs11_21 =  zvis11 * akappa(ji,jj-1,1,2) + zdiag 
     341               zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1) 
     342               zs22_21 =  zvis11 * akappa(ji,jj-1,2,2) + zdiag 
     343               zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1) 
     344 
     345               zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
     346               zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
     347               zvis12 =       zviseta (ji-1,jj) + dm 
     348               zvis21 =       zviseta (ji-1,jj) 
     349               zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) ) 
     350               zs11_12 =  zvis11 * akappa(ji-1,jj,1,2) + zdiag 
     351               zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1) 
     352               zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag 
     353               zs21_12 =  zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1) 
     354 
     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               zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) ) 
     360               zs11_22 =  zvis11 * akappa(ji,jj,1,2) + zdiag 
     361               zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1) 
     362               zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag 
     363               zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1) 
     364 
     365               zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22   & 
     366                  &          - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12   & 
     367                  &          - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11   & 
     368                  &          + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12   & 
     369                  &          + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21   & 
     370                  &          + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22   & 
     371                  &          - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21   & 
     372                  &          - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 
     373 
     374               zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22   & 
     375                  &          - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12   & 
     376                  &          - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11   & 
     377                  &          + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12   & 
     378                  &          - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21   & 
     379                  &          - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22   & 
     380                  &          + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21   & 
     381                  &          + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 
    382382            END DO 
    383  
     383         END DO 
     384 
     385         ! GAUSS-SEIDEL method 
     386         !                                                      ! ================ ! 
     387iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! 
     388            !                                                   ! ================ ! 
     389!CDIR NOVERRCHK 
    384390            DO jj = k_j1+1, k_jpj-1 
    385                zsang  = SIGN( 1.e0, fcor(1,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere 
    386                DO ji = 2, jpim1 
    387                  zur     = u_ice(ji,jj) - u_io(ji,jj) 
    388                  zvr     = v_ice(ji,jj) - v_io(ji,jj) 
    389                  zmod    = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) ) 
    390                  za      = rhoco * zmod 
    391                  zac     = za * cangvg 
    392                   zmpzas  = alpha * zcorl(ji,jj) + za * zsang 
     391!CDIR NOVERRCHK 
     392               DO ji = fs_2, fs_jpim1 
     393                  ! 
     394                  ze11 =   akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj) 
     395                  ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj) 
     396                  ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj) 
     397                  ze21 =   akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj) 
     398                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
     399                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
     400                  zvis12 =       zviseta (ji,jj-1) + dm 
     401                  zvis21 =       zviseta (ji,jj-1) 
     402                  zdiag = zvis22 * ( ze11 + ze22 ) 
     403                  zs11_21 =  zvis11 * ze11 + zdiag 
     404                  zs12_21 =  zvis12 * ze12 + zvis21 * ze21 
     405                  zs22_21 =  zvis11 * ze22 + zdiag 
     406                  zs21_21 =  zvis12 * ze21 + zvis21 * ze12 
     407 
     408                  ze11 =   akappa(ji-1,jj,1,1) * ( zu_a(ji  ,jj+1) - zu_a(ji-1,jj+1) )   & 
     409                     &   + akappa(ji-1,jj,1,2) * ( zv_a(ji  ,jj+1) + zv_a(ji-1,jj+1) ) 
     410                  ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji  ,jj+1) )   & 
     411                     &   - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji  ,jj+1) ) 
     412                  ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji  ,jj+1) )   & 
     413                     &   + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji  ,jj+1) ) 
     414                  ze21 =   akappa(ji-1,jj,1,1) * ( zv_a(ji  ,jj+1) - zv_a(ji-1,jj+1) )   & 
     415                     &   - akappa(ji-1,jj,1,2) * ( zu_a(ji  ,jj+1) + zu_a(ji-1,jj+1) ) 
     416                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
     417                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
     418                  zvis12 =       zviseta (ji-1,jj) + dm 
     419                  zvis21 =       zviseta (ji-1,jj) 
     420                  zdiag = zvis22 * ( ze11 + ze22 ) 
     421                  zs11_12 =  zvis11 * ze11 + zdiag 
     422                  zs12_12 =  zvis12 * ze12 + zvis21 * ze21 
     423                  zs22_12 =  zvis11 * ze22 + zdiag 
     424                  zs21_12 =  zvis12 * ze21 + zvis21 * ze12 
     425 
     426                  ze11 =   akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji  ,jj+1) )   & 
     427                     &   + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji  ,jj+1) ) 
     428                  ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji  ,jj+1) - zu_a(ji+1,jj+1) )   & 
     429                     &   - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji  ,jj+1) + zv_a(ji+1,jj+1) ) 
     430                  ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji  ,jj+1) - zv_a(ji+1,jj+1) )   & 
     431                     &   + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji  ,jj+1) + zu_a(ji+1,jj+1) ) 
     432                  ze21 =   akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji  ,jj+1) )   & 
     433                     &   - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji  ,jj+1) ) 
     434                  zvis11 = 2.0 * zviseta (ji,jj) + dm 
     435                  zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj) 
     436                  zvis12 =       zviseta (ji,jj) + dm 
     437                  zvis21 =       zviseta (ji,jj) 
     438                  zdiag = zvis22 * ( ze11 + ze22 ) 
     439                  zs11_22 =  zvis11 * ze11 + zdiag 
     440                  zs12_22 =  zvis12 * ze12 + zvis21 * ze21 
     441                  zs22_22 =  zvis11 * ze22 + zdiag 
     442                  zs21_22 =  zvis12 * ze21 + zvis21 * ze12 
     443 
     444            ! 2nd part 
     445                  ze11 =   akappa(ji-1,jj-1,1,1) * ( zu_a(ji  ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) )   & 
     446                     &   + akappa(ji-1,jj-1,1,2) * ( zv_a(ji  ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) ) 
     447                  ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji  ,jj-1) - zu_a(ji-1,jj) )   & 
     448                     &   - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji  ,jj-1) + zv_a(ji-1,jj) ) 
     449                  ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji  ,jj-1) - zv_a(ji-1,jj) )   & 
     450                     &   + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji  ,jj-1) + zu_a(ji-1,jj) ) 
     451                  ze21 =   akappa(ji-1,jj-1,1,1) * ( zv_a(ji  ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) )   & 
     452                     &  -  akappa(ji-1,jj-1,1,2) * ( zu_a(ji  ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) ) 
     453                  zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 
     454                  zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 
     455                  zvis12 =       zviseta (ji-1,jj-1) + dm 
     456                  zvis21 =       zviseta (ji-1,jj-1) 
     457                  zdiag = zvis22 * ( ze11 + ze22 ) 
     458                  zs11_11 =  zvis11 * ze11 + zdiag 
     459                  zs12_11 =  zvis12 * ze12 + zvis21 * ze21 
     460                  zs22_11 =  zvis11 * ze22 + zdiag 
     461                  zs21_11 =  zvis12 * ze21 + zvis21 * ze12 
     462 
     463                  ze11 =   akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji  ,jj-1) )   & 
     464                     &   + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji  ,jj-1) ) 
     465                  ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji  ,jj-1) + zu_a(ji+1,jj-1) )   & 
     466                     &   - akappa(ji,jj-1,2,1) * ( zv_a(ji  ,jj-1) + zv_a(ji+1,jj-1) ) 
     467                  ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji  ,jj-1) + zv_a(ji+1,jj-1) )   & 
     468                     &   + akappa(ji,jj-1,2,1) * ( zu_a(ji  ,jj-1) + zu_a(ji+1,jj-1) ) 
     469                  ze21 =   akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji  ,jj-1) )   & 
     470                     &   - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji  ,jj-1) ) 
     471                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm 
     472                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1) 
     473                  zvis12 =       zviseta (ji,jj-1) + dm 
     474                  zvis21 =       zviseta (ji,jj-1) 
     475                  zdiag = zvis22 * ( ze11 + ze22 ) 
     476                  zs11_21 =  zs11_21 + zvis11 * ze11 + zdiag 
     477                  zs12_21 =  zs12_21 + zvis12 * ze12 + zvis21 * ze21 
     478                  zs22_21 =  zs22_21 + zvis11 * ze22 + zdiag 
     479                  zs21_21 =  zs21_21 + zvis12 * ze21 + zvis21 * ze12 
     480 
     481                  ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj) 
     482                  ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj) 
     483                  ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj) 
     484                  ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj) 
     485                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm 
     486                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj) 
     487                  zvis12 =       zviseta (ji-1,jj) + dm 
     488                  zvis21 =       zviseta (ji-1,jj) 
     489                  zdiag = zvis22 * ( ze11 + ze22 ) 
     490                  zs11_12 =  zs11_12 + zvis11 * ze11 + zdiag 
     491                  zs12_12 =  zs12_12 + zvis12 * ze12 + zvis21 * ze21 
     492                  zs22_12 =  zs22_12 + zvis11 * ze22 + zdiag 
     493                  zs21_12 =  zs21_12 + zvis12 * ze21 + zvis21 * ze12 
     494 
     495                  zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22  & 
     496                     &  - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12  & 
     497                     &  - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11  & 
     498                     &  + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12  & 
     499                     &  + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21  & 
     500                     &  + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22  & 
     501                     &  - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21  & 
     502                     &  - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 
     503 
     504                  zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22  & 
     505                     &  - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12  & 
     506                     &  - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11  & 
     507                     &  + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12  & 
     508                     &  - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21  & 
     509                     &  - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22  & 
     510                     &  + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21  & 
     511                     &  + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 
     512 
     513                  zur     = zu_a(ji,jj) - ui_oce(ji,jj) 
     514                  zvr     = zv_a(ji,jj) - vi_oce(ji,jj) 
     515!!!! 
     516                  zmod    = SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 
     517                  za      = rhoco * zmod 
     518!!!! 
     519!!gm chg resul    za      = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 
     520                  zac     = za * cangvg 
     521                  zmpzas  = alpha * zcorl(ji,jj) + za * zsang(ji,jj) 
    393522                  zmassdt = zusdtp * zmass(ji,jj) 
    394523                  zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj) 
    395524 
    396                   za1(ji,jj) =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   & 
    397                      &        + za * ( cangvg * u_io(ji,jj) - zsang * v_io(ji,jj) ) 
    398  
    399                   za2(ji,jj) =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   & 
    400                      &        + za * ( cangvg * v_io(ji,jj) + zsang * u_io(ji,jj) ) 
    401  
    402                   zb1(ji,jj)  = zmassdt + zac - zc1u(ji,jj) 
    403                   zb2(ji,jj)  = zmpzas  - zc2u(ji,jj) 
    404                   zc1(ji,jj)  = zmpzas  + zc1v(ji,jj) 
    405                   zc2(ji,jj)  = zmassdt + zac  - zc2v(ji,jj)  
    406                   zdeter      = zc1(ji,jj) * zb2(ji,jj) + zc2(ji,jj) * zb1(ji,jj) 
    407                   zden(ji,jj) = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 
     525                  za1 =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   & 
     526                     &        + za * ( cangvg * ui_oce(ji,jj) - zsang(ji,jj) * vi_oce(ji,jj) ) 
     527                  za2 =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   & 
     528                     &        + za * ( cangvg * vi_oce(ji,jj) + zsang(ji,jj) * ui_oce(ji,jj) ) 
     529                  zb1    = zmassdt + zac - zc1u(ji,jj) 
     530                  zb2    = zmpzas        - zc2u(ji,jj) 
     531                  zc1    = zmpzas        + zc1v(ji,jj) 
     532                  zc2    = zmassdt + zac - zc2v(ji,jj) 
     533                  zdeter = zc1 * zb2 + zc2 * zb1 
     534                  zden   = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 
     535                  zunw   = (  ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden 
     536                  zvnw   = (  ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden 
     537                  zmask  = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 
     538 
     539                  zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask 
     540                  zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask 
    408541               END DO 
    409542            END DO 
    410543 
    411             ! The computation of ice interaction term is splitted into two parts 
    412             !------------------------------------------------------------------------- 
    413  
    414             ! Terms that do not involve already up-dated velocities. 
    415           
    416             DO jj = k_j1+1, k_jpj-1 
    417                DO ji = 2, jpim1 
    418                   iim1 = ji 
    419                   ijm1 = jj - 1 
    420                   iip1 = ji + 1 
    421                   ijp1 = jj 
    422                   ze11 =   akappa(iim1,ijm1,1,1) * u_ice(iip1,ijp1) + akappa(iim1,ijm1,1,2) * v_ice(iip1,ijp1) 
    423                   ze12 = + akappa(iim1,ijm1,2,2) * u_ice(iip1,ijp1) - akappa(iim1,ijm1,2,1) * v_ice(iip1,ijp1) 
    424                   ze22 = + akappa(iim1,ijm1,2,2) * v_ice(iip1,ijp1) + akappa(iim1,ijm1,2,1) * u_ice(iip1,ijp1) 
    425                   ze21 =   akappa(iim1,ijm1,1,1) * v_ice(iip1,ijp1) - akappa(iim1,ijm1,1,2) * u_ice(iip1,ijp1) 
    426                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    427                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    428                   zvis12 =       zviseta (iim1,ijm1) + dm 
    429                   zvis21 =       zviseta (iim1,ijm1) 
    430                   zdiag = zvis22 * ( ze11 + ze22 ) 
    431                   zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag 
    432                   zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21 
    433                   zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag 
    434                   zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12 
    435  
    436  
    437                   iim1 = ji - 1 
    438                   ijm1 = jj 
    439                   iip1 = ji 
    440                   ijp1 = jj + 1                    
    441                   ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   & 
    442                      &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijp1) + v_ice(iim1,ijp1) ) 
    443                   ze12 = + akappa(iim1,ijm1,2,2) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )   & 
    444                      &   - akappa(iim1,ijm1,2,1) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) 
    445                   ze22 = + akappa(iim1,ijm1,2,2) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )   & 
    446                      &   + akappa(iim1,ijm1,2,1) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) 
    447                   ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   & 
    448                      &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 
    449                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    450                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    451                   zvis12 =       zviseta (iim1,ijm1) + dm 
    452                   zvis21 =       zviseta (iim1,ijm1) 
    453                   zdiag = zvis22 * ( ze11 + ze22 ) 
    454                   zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag 
    455                   zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21 
    456                   zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag 
    457                   zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12 
    458  
    459                   iim1 = ji 
    460                   ijm1 = jj 
    461                   iip1 = ji + 1 
    462                   ijp1 = jj + 1 
    463                   ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   & 
    464                      &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) + v_ice(iim1,ijp1) ) 
    465                   ze12 = - akappa(iim1,ijm1,2,2) * ( u_ice(iip1,ijm1) - u_ice(iim1,ijp1) - u_ice(iip1,ijp1) )   & 
    466                      &   - akappa(iim1,ijm1,2,1) * ( v_ice(iip1,ijm1) + v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) 
    467                   ze22 = - akappa(iim1,ijm1,2,2) * ( v_ice(iip1,ijm1) - v_ice(iim1,ijp1) - v_ice(iip1,ijp1) )   & 
    468                      &   + akappa(iim1,ijm1,2,1) * ( u_ice(iip1,ijm1) + u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) 
    469                   ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   & 
    470                      &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) + u_ice(iim1,ijp1) )  
    471                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    472                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    473                   zvis12 =       zviseta (iim1,ijm1) + dm 
    474                   zvis21 =       zviseta (iim1,ijm1) 
    475  
    476                   zdiag = zvis22 * ( ze11 + ze22 ) 
    477                   zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag 
    478                   zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21 
    479                   zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag 
    480                   zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12 
    481  
    482                END DO 
     544            CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. ) 
     545            CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. ) 
     546 
     547            ! Test of Convergence 
     548            DO jj = k_j1+1 , k_jpj-1 
     549               zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) ) 
    483550            END DO 
    484  
    485             ! Terms involving already up-dated velocities. 
    486             !-Using the arrays zu_ice and zv_ice in the computation of the terms ze leads to JACOBI's method;  
    487             ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method. 
    488               
    489             DO jj = k_j1+1, k_jpj-1 
    490                DO ji = 2, jpim1 
    491                   iim1 = ji - 1 
    492                   ijm1 = jj - 1 
    493                   iip1 = ji 
    494                   ijp1 = jj 
    495                   ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) - zu_ice(iim1,ijp1) )   & 
    496                      &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) + zv_ice(iim1,ijp1) ) 
    497                   ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) - zu_ice(iim1,ijp1) )   & 
    498                      &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) + zv_ice(iim1,ijp1) ) 
    499                   ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) - zv_ice(iim1,ijp1) )   & 
    500                      &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) + zu_ice(iim1,ijp1) ) 
    501                   ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) - zv_ice(iim1,ijp1) )   & 
    502                      &  -  akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) + zu_ice(iim1,ijp1) ) 
    503                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    504                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    505                   zvis12 =       zviseta (iim1,ijm1) + dm 
    506                   zvis21 =       zviseta (iim1,ijm1) 
    507  
    508                   zdiag = zvis22 * ( ze11 + ze22 ) 
    509                   zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag 
    510                   zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21 
    511                   zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag 
    512                   zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12 
    513  
    514 #if defined key_agrif 
    515              END DO 
    516           END DO 
    517  
    518           DO jj = k_j1+1, k_jpj-1 
    519              DO ji = 2, jpim1 
    520 #endif 
    521  
    522                   iim1 = ji 
    523                   ijm1 = jj - 1 
    524                   iip1 = ji + 1 
    525                   ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) )   & 
    526                      &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) ) 
    527                   ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )   & 
    528                      &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) ) 
    529                   ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )   & 
    530                      &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) ) 
    531                   ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) )   & 
    532                      &   - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) ) 
    533                   zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 
    534                   zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 
    535                   zvis12 =       zviseta (iim1,ijm1) + dm 
    536                   zvis21 =       zviseta (iim1,ijm1) 
    537  
    538                   zdiag = zvis22 * ( ze11 + ze22 ) 
    539                   zs11(ji,jj,2,1) =  zs11(ji,jj,2,1) + zvis11 * ze11 + zdiag 
    540                   zs12(ji,jj,2,1) =  zs12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21 
    541                   zs22(ji,jj,2,1) =  zs22(ji,jj,2,1) + zvis11 * ze22 + zdiag 
    542                   zs21(ji,jj,2,1) =  zs21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12 
    543  
    544  
    545                   iim1 = ji - 1 
    546                   ijm1 = jj  
    547                   ze11 = - akappa(iim1,ijm1,1,1) * zu_ice(iim1,ijm1) + akappa(iim1,ijm1,1,2) * zv_ice(iim1,ijm1) 
    548                   ze12 = - akappa(iim1,ijm1,2,2) * zu_ice(iim1,ijm1) - akappa(iim1,ijm1,2,1) * zv_ice(iim1,ijm1) 
    549                   ze22 = - akappa(iim1,ijm1,2,2) * zv_ice(iim1,ijm1) + akappa(iim1,ijm1,2,1) * zu_ice(iim1,ijm1) 
    550                   ze21 = - akappa(iim1,ijm1,1,1) * zv_ice(iim1,ijm1) - akappa(iim1,ijm1,1,2) * 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,1,2) =  zs11(ji,jj,1,2) + zvis11 * ze11 + zdiag  
    558                   zs12(ji,jj,1,2) =  zs12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21 
    559                   zs22(ji,jj,1,2) =  zs22(ji,jj,1,2) + zvis11 * ze22 + zdiag 
    560                   zs21(ji,jj,1,2) =  zs21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12 
    561  
    562 #if defined key_agrif 
    563              END DO 
    564           END DO 
    565  
    566           DO jj = k_j1+1, k_jpj-1 
    567              DO ji = 2, jpim1 
    568 #endif 
    569                   zd1(ji,jj) =   & 
    570                      + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)  & 
    571                      - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)  & 
    572                      - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)  & 
    573                      + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)  & 
    574                      + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)  & 
    575                      + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)  & 
    576                      - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)  & 
    577                      - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 
    578                   zd2(ji,jj) =   & 
    579                      + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)  & 
    580                      - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)  & 
    581                      - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)  & 
    582                      + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)  & 
    583                      - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)  & 
    584                      - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)  & 
    585                      + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)  & 
    586                      + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 
    587                END DO 
     551            zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) ) 
     552!!!! this should be faster on scalar processor 
     553!           zresm = MAXVAL(  MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ),   & 
     554!              &                  ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) )  ) 
     555!!!! 
     556            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
     557 
     558            DO jj = k_j1, k_jpj 
     559               zu_a(:,jj) = zu_n(:,jj) 
     560               zv_a(:,jj) = zv_n(:,jj) 
    588561            END DO 
    589562 
    590             DO jj = k_j1+1, k_jpj-1 
    591                DO ji = 2, jpim1 
    592                   zunw = (  ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj)        & 
    593                      &    + ( za2(ji,jj) + zd2(ji,jj) ) * zc1(ji,jj) ) * zden(ji,jj) 
    594  
    595                   zvnw = (  ( za2(ji,jj) + zd2(ji,jj) ) * zb1(ji,jj)        & 
    596                      &    - ( za1(ji,jj) + zd1(ji,jj) ) * zb2(ji,jj) ) * zden(ji,jj) 
    597  
    598                   zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 
    599  
    600                   u_ice(ji,jj) = ( u_ice(ji,jj) + om * ( zunw - u_ice(ji,jj) ) * tmu(ji,jj) ) * zmask 
    601                   v_ice(ji,jj) = ( v_ice(ji,jj) + om * ( zvnw - v_ice(ji,jj) ) * tmu(ji,jj) ) * zmask 
    602                END DO 
     563            IF( zresm <= resl )   EXIT   iflag 
     564 
     565            !                                                   ! ================ ! 
     566         END DO    iflag                                        !  end Relaxation  ! 
     567         !                                                      ! ================ ! 
     568 
     569         IF( zindu == 0 ) THEN      ! even iteration 
     570            DO jj = k_j1 , k_jpj-1 
     571               zu0(:,jj) = zu_a(:,jj) 
     572               zv0(:,jj) = zv_a(:,jj) 
    603573            END DO 
    604  
    605             CALL lbc_lnk( u_ice, 'I', -1. ) 
    606             CALL lbc_lnk( v_ice, 'I', -1. ) 
    607  
    608             !---  5.2.5.4. Convergence test. 
    609             DO jj = k_j1+1 , k_jpj-1 
    610                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    611             END DO 
    612             zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 
    613             IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
    614  
    615             IF ( zresm <= resl) EXIT iflag 
    616  
    617          END DO iflag 
    618  
    619          zindu1 = 1.0 - zindu 
    620          DO jj = k_j1 , k_jpj-1 
    621             zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj) 
    622             zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj) 
    623          END DO 
    624       !                                                   ! ==================== ! 
     574         ENDIF 
     575         !                                                ! ==================== ! 
    625576      END DO                                              !  end loop over iter  ! 
    626577      !                                                   ! ==================== ! 
     578 
     579      ui_ice(:,:) = zu_a(:,1:jpj) 
     580      vi_ice(:,:) = zv_a(:,1:jpj) 
    627581 
    628582      IF(ln_ctl) THEN 
    629583         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter 
    630584         CALL prt_ctl_info(charout) 
    631          CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 
     585         CALL prt_ctl(tab2d_1=ui_ice, clinfo1=' lim_rhg  : ui_ice :', tab2d_2=vi_ice, clinfo2=' vi_ice :') 
    632586      ENDIF 
    633587 

  • trunk/NEMO/LIM_SRC/limrst.F90

    r699 r717  
    1717   !!---------------------------------------------------------------------- 
    1818   USE ice 
    19    USE dom_oce 
    20    USE ice_oce         ! ice variables 
     19   USE sbc_oce 
     20   USE sbc_ice 
    2121   USE daymod 
    2222 
    … …  
    2828 
    2929   PUBLIC   lim_rst_opn     ! routine called by ??? module 
    30    PUBLIC   lim_rst_write   ! routine called by ??? module 
    31    PUBLIC   lim_rst_read    ! routine called by ??? module 
     30   PUBLIC   lim_rst_write   ! routine called by lim_step.F90 
     31   PUBLIC   lim_rst_read    ! routine called by lim_init.F90 
    3232 
    3333   LOGICAL, PUBLIC ::   lrst_ice         !: logical to control the ice restart write  
    … …  
    5656      IF( kt == nit000 )   lrst_ice = .FALSE. 
    5757       
    58       IF( kt == nitrst - 2*nfice + 1 .OR.  nitend - nit000 + 1 <= nfice ) THEN 
    59          ! beware if model runs less than nfice + 1 time step 
     58      IF( kt == nitrst - 2*nn_fsbc + 1 .OR.  nitend - nit000 + 1 <= nn_fsbc ) THEN 
     59         ! beware if model runs less than nn_fsbc + 1 time step 
    6060         ! beware of the format used to write kt (default is i8.8, that should be large enough) 
    6161         IF( nitrst > 1.0e9 ) THEN    
    … …  
    7474   END SUBROUTINE lim_rst_opn 
    7575 
     76 
    7677   SUBROUTINE lim_rst_write( kt ) 
    7778      !!---------------------------------------------------------------------- 
    … …  
    8081      !! ** purpose  :   output of sea-ice variable in a netcdf file 
    8182      !!---------------------------------------------------------------------- 
    82       INTEGER, INTENT(in) ::   kt     ! number of iteration 
    83       !! 
    84       INTEGER ::   iter     ! kt + nfice -1 
    85       !!---------------------------------------------------------------------- 
    86  
    87       iter = kt + nfice -1 
     83      INTEGER, INTENT(in) ::   kt   ! number of iteration 
     84      ! 
     85      INTEGER ::   iter   ! kt + nn_fsbc -1 
     86      !!---------------------------------------------------------------------- 
     87 
     88      iter = kt + nn_fsbc - 1 
    8889 
    8990      IF( iter == nitrst ) THEN 
    9091         IF(lwp) WRITE(numout,*) 
    9192         IF(lwp) WRITE(numout,*) 'lim_rst_write : write ice restart.output NetCDF file  kt =', kt 
    92          IF(lwp) WRITE(numout,*) '~~~~~~~'          
    93       ENDIF 
    94  
    95       ! Write in numriw (if iter == nitrst) 
     93         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 
     94      ENDIF 
     95 
     96      ! Write in numriw (if iter  == nitrst) 
    9697      ! ------------------  
    9798      !                                                                     ! calendar control 
    98       CALL iom_rstput( iter, nitrst, numriw, 'nfice' , REAL( nfice, wp) )      ! time-step  
    99       CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) )      ! date 
     99      CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) )  
    100100       
    101101      CALL iom_rstput( iter, nitrst, numriw, 'hicif' , hicif (:,:)   )      ! prognostic variables  
    … …  
    109109      CALL iom_rstput( iter, nitrst, numriw, 'tbif2' , tbif  (:,:,2) ) 
    110110      CALL iom_rstput( iter, nitrst, numriw, 'tbif3' , tbif  (:,:,3) ) 
    111       CALL iom_rstput( iter, nitrst, numriw, 'u_ice' , u_ice (:,:)   ) 
    112       CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice (:,:)   ) 
    113       CALL iom_rstput( iter, nitrst, numriw, 'gtaux' , gtaux (:,:)   ) 
    114       CALL iom_rstput( iter, nitrst, numriw, 'gtauy' , gtauy (:,:)   ) 
     111      CALL iom_rstput( iter, nitrst, numriw, 'ui_ice', ui_ice(:,:)   ) 
     112      CALL iom_rstput( iter, nitrst, numriw, 'vi_ice', vi_ice(:,:)   ) 
    115113      CALL iom_rstput( iter, nitrst, numriw, 'qstoif', qstoif(:,:)   ) 
    116114      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq (:,:)   ) 
    … …  
    165163      !! ** purpose  :   read of sea-ice variable restart in a netcdf file 
    166164      !!---------------------------------------------------------------------- 
    167       ! 
    168       REAL(wp) ::   zfice, ziter 
     165      REAL(wp) ::   ziter 
    169166      !!---------------------------------------------------------------------- 
    170167 
    … …  
    172169         WRITE(numout,*) 
    173170         WRITE(numout,*) 'lim_rst_read : read ice NetCDF restart file' 
    174          WRITE(numout,*) '~~~~~~~~' 
     171         WRITE(numout,*) '~~~~~~~~~~~~' 
    175172      ENDIF 
    176173 
    177174      CALL iom_open ( 'restart_ice_in', numrir, kiolib = jprstlib ) 
    178175 
    179       CALL iom_get( numrir, 'nfice' , zfice ) 
    180       CALL iom_get( numrir, 'kt_ice', ziter )     
    181       IF(lwp) WRITE(numout,*) '   read ice restart file at time step    : ', ziter 
     176      CALL iom_get( numrir, 'kt_ice' , ziter ) 
     177      IF(lwp) WRITE(numout,*) '   read ice restart file at time step    : ', INT( ziter ) 
    182178      IF(lwp) WRITE(numout,*) '   in any case we force it to nit000 - 1 : ', nit000 - 1 
    183179 
    … …  
    186182      IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 )   & 
    187183         &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart',  & 
    188          &                   '   verify the file or rerun with the value 0 for the',        & 
    189          &                   '   control of time parameter  nrstdt' ) 
    190       IF( INT(zfice) /= nfice          .AND. ABS( nrstdt ) == 1 )   & 
    191          &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nfice in ice restart',  & 
    192184         &                   '   verify the file or rerun with the value 0 for the',        & 
    193185         &                   '   control of time parameter  nrstdt' ) 
    … …  
    203195      CALL iom_get( numrir, jpdom_autoglo, 'tbif2' , tbif(:,:,2) )     
    204196      CALL iom_get( numrir, jpdom_autoglo, 'tbif3' , tbif(:,:,3) )     
    205       CALL iom_get( numrir, jpdom_autoglo, 'u_ice' , u_ice  )     
    206       CALL iom_get( numrir, jpdom_autoglo, 'v_ice' , v_ice  )     
    207       CALL iom_get( numrir, jpdom_autoglo, 'gtaux' , gtaux  )     
    208       CALL iom_get( numrir, jpdom_autoglo, 'gtauy' , gtauy  )     
     197      CALL iom_get( numrir, jpdom_autoglo, 'ui_ice', ui_ice )     
     198      CALL iom_get( numrir, jpdom_autoglo, 'vi_ice', vi_ice )     
    209199      CALL iom_get( numrir, jpdom_autoglo, 'qstoif', qstoif )     
    210200      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq  )     

  • trunk/NEMO/LIM_SRC/limsbc.F90

    r702 r717  
    6161      !!              - Update  
    6262      !!      
    63       !! ** Outputs : - fsolar  : solar heat flux at sea ice/ocean interface 
    64       !!              - fnsolar : non solar heat flux  
    65       !!              - fsalt   : salt flux at sea ice/ocean interface 
    66       !!              - fmass   : freshwater flux at sea ice/ocean interface 
     63      !! ** Outputs : - qsr    : sea heat flux:     solar  
     64      !!              - qns    : sea heat flux: non solar 
     65      !!              - emp    : freshwater budget: volume flux  
     66      !!              - emps   : freshwater budget: concentration/dillution  
     67      !!              - utau   : sea surface i-stress (ocean referential) 
     68      !!              - vtau   : sea surface j-stress (ocean referential) 
    6769      !! 
    6870      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 
    … …  
    183185            DO ji = 1, jpi 
    184186               ! ... change the cosinus angle sign in the south hemisphere 
    185                zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg  
     187!!              zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg    ! do the full loop and avoid lbc_lnk 
     188               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg 
    186189               ! ... ice velocity relative to the ocean 
    187                zu_io  = u_ice(ji,jj) - u_oce(ji,jj) 
    188                zv_io  = v_ice(ji,jj) - v_oce(ji,jj) 
    189                zmod   = SQRT( zu_io * zu_io + zv_io * zv_io ) 
     190               zu_io  = ui_ice(ji,jj) - ui_oce(ji,jj) 
     191               zv_io  = vi_ice(ji,jj) - vi_oce(ji,jj) 
     192               zmod   = rhoco * SQRT( zu_io * zu_io + zv_io * zv_io ) 
    190193               ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point) 
    191                ztio_u(ji,jj) = cw * zmod * ( cangvg * zu_io - zsang * zv_io ) 
    192                ztio_v(ji,jj) = cw * zmod * ( cangvg * zv_io + zsang * zu_io ) 
     194               ztio_u(ji,jj) = zmod * ( cangvg * zu_io - zsang * zv_io ) 
     195               ztio_v(ji,jj) = zmod * ( cangvg * zv_io + zsang * zu_io ) 
    193196               !  
    194197            END DO 
    … …  
    196199 
    197200         DO jj = 2, jpjm1 
    198             DO ji = 2, jpim1 
     201            DO ji = fs_2, fs_jpim1   ! vertor opt. 
    199202               ! ... ice-cover wheighted ice-ocean stress at U and V-points  (from I-point values) 
    200203               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) ) 

  • trunk/NEMO/LIM_SRC/limthd.F90

    r709 r717  
    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_ice_lim 
    711   !!---------------------------------------------------------------------- 
    … …  
    1115   !!   lim_thd_init : initialisation of sea-ice thermodynamic 
    1216   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
     17   USE phycst          ! physical constants 
    1418   USE dom_oce         ! ocean space and time domain variables 
    15    USE sbc_oce         ! surface boundary condition: ocean 
     19   USE lbclnk 
     20   USE in_out_manager  ! I/O manager 
     21   USE ice             ! LIM sea-ice variables 
    1622   USE ice_oce         ! sea-ice/ocean variables 
     23   USE sbc_oce         !  
     24   USE sbc_ice         !  
    1725   USE thd_ice         ! LIM thermodynamic sea-ice variables 
    1826   USE dom_ice         ! LIM sea-ice domain 
    19    USE ice             ! LIM sea-ice variables 
    2027   USE iceini 
    2128   USE limthd_zdf 
    2229   USE limthd_lac 
    2330   USE limtab 
    24    USE phycst          ! physical constants 
    25    USE in_out_manager  ! I/O manager 
    2631   USE prtctl          ! Print control 
    27    USE lbclnk 
    2832       
    2933   IMPLICIT NONE 
    3034   PRIVATE 
    3135 
    32    !! * Routine accessibility 
    33    PUBLIC lim_thd       ! called by lim_step 
    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 
     36   PUBLIC   lim_thd    ! called by lim_step 
     37 
     38   REAL(wp)  ::   epsi20 = 1.e-20   ,  &  ! constant values 
     39      &           epsi16 = 1.e-16   ,  & 
     40      &           epsi04 = 1.e-04   ,  & 
     41      &           zzero  = 0.e0     ,  & 
     42      &           zone   = 1.e0 
    4243 
    4344   !! * Substitutions 
    … …  
    4748   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    4849   !! $Id$ 
    49    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     50   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5051   !!---------------------------------------------------------------------- 
    5152 
    … …  
    6869      !!             - back to the geographic grid 
    6970      !! 
    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 
     71      !! References :   Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 
    7672      !!--------------------------------------------------------------------- 
    7773      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    78  
     74      !! 
    7975      INTEGER  ::   ji, jj,    &   ! dummy loop indices 
    8076         nbpb  ,               &   ! nb of icy pts for thermo. cal. 
    … …  
    9288         zfontn             ,  &   ! heat flux from snow thickness 
    9389         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 
     90      REAL(wp), DIMENSION(jpi,jpj) ::   zhicifp,   &  ! ice thickness for outputs 
     91         &                              zqlbsbq       ! link with lead energy budget qldif 
     92      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmsk      ! working array 
    9993      !!------------------------------------------------------------------- 
    10094 
    101       IF( kt == nit000  )   CALL lim_thd_init  ! Initialization (first time-step only) 
     95      IF( kt == nit000 )   CALL lim_thd_init  ! Initialization (first time-step only) 
    10296    
    10397      !-------------------------------------------! 
    … …  
    188182            !  temperature and turbulent mixing (McPhee, 1992) 
    189183            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) )  
     184!!gm old    fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_io(ji,jj) - tfu(ji,jj) )  
     185            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) )  
    191186            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice 
    192187                         
    193188            !  partial computation of the lead energy budget (qldif) 
    194189            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) )   & 
     190            zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean 
     191            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   & 
    197192               &                               + zfnsol + fdtcn(ji,jj) - zfontn     & 
    198193               &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   & 
    … …  
    206201             
    207202            !  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 ) 
     203!!gm old    qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_io(ji,jj) ) * ( 1 - zinda ) 
     204            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda ) 
    209205             
    210206            !  calculate oceanic heat flux. 
    … …  
    258254         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb)     , fr1_i0     , jpi, jpj, npb(1:nbpb) ) 
    259255         CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb)     , fr2_i0     , jpi, jpj, npb(1:nbpb) ) 
    260          CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb)     , qnsr_ice   , jpi, jpj, npb(1:nbpb) ) 
     256         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb)     , qns_ice    , jpi, jpj, npb(1:nbpb) ) 
    261257#if ! defined key_coupled 
    262258         CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb)     , qla_ice    , jpi, jpj, npb(1:nbpb) ) 
    … …  
    394390      IF(ln_ctl) THEN 
    395391         CALL prt_ctl_info(' lim_thd  end  ') 
    396          CALL prt_ctl(tab2d_1=hicif , clinfo1=' lim_thd: hicif   : ', tab2d_2=hsnif , clinfo2=' hsnif  : ') 
    397          CALL prt_ctl(tab2d_1=frld  , clinfo1=' lim_thd: frld    : ', tab2d_2=hicifp, clinfo2=' hicifp : ') 
    398          CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif  : ', tab2d_2=pfrld , clinfo2=' pfrld  : ') 
    399          CALL prt_ctl(tab2d_1=sist  , clinfo1=' lim_thd: sist    : ') 
    400          CALL prt_ctl(tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1  : ') 
    401          CALL prt_ctl(tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2  : ') 
    402          CALL prt_ctl(tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3  : ') 
    403          CALL prt_ctl(tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn   : ', tab2d_2=qdtcn , clinfo2=' qdtcn  : ') 
    404          CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif  : ', tab2d_2=fsbbq , clinfo2=' fsbbq  : ') 
    405       ENDIF 
    406  
    407     END SUBROUTINE lim_thd 
    408  
    409  
    410     SUBROUTINE lim_thd_init 
     392         CALL prt_ctl(tab2d_1=hicif      , clinfo1=' hicif   : ', tab2d_2=hsnif      , clinfo2=' hsnif  : ') 
     393         CALL prt_ctl(tab2d_1=frld       , clinfo1=' frld    : ', tab2d_2=hicifp     , clinfo2=' hicifp : ') 
     394         CALL prt_ctl(tab2d_1=phicif     , clinfo1=' phicif  : ', tab2d_2=pfrld      , clinfo2=' pfrld  : ') 
     395         CALL prt_ctl(tab2d_1=sist       , clinfo1=' sist    : ', tab2d_2=tbif(:,:,1), clinfo2=' tbif 1 : ') 
     396         CALL prt_ctl(tab2d_1=tbif(:,:,2), clinfo1=' tbif 2  : ', tab2d_2=tbif(:,:,3), clinfo2=' tbif 3 : ') 
     397         CALL prt_ctl(tab2d_1=fdtcn      , clinfo1=' fdtcn   : ', tab2d_2=qdtcn      , clinfo2=' qdtcn  : ') 
     398         CALL prt_ctl(tab2d_1=qstoif     , clinfo1=' qstoif  : ', tab2d_2=fsbbq      , clinfo2=' fsbbq  : ') 
     399      ENDIF 
     400       ! 
     401   END SUBROUTINE lim_thd 
     402 
     403 
     404   SUBROUTINE lim_thd_init 
    411405      !!------------------------------------------------------------------- 
    412406      !!                   ***  ROUTINE lim_thd_init ***  
    … …  
    419413      !! 
    420414      !! ** input   :   Namelist namicether 
    421       !! 
    422       !! history : 
    423       !!  8.5  ! 03-08 (C. Ethe) original code 
    424415      !!------------------------------------------------------------------- 
    425416      NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax  ,        & 

  • trunk/NEMO/LIM_SRC/limthd_zdf.F90

    r699 r717  
    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_ice_lim 
    710   !!---------------------------------------------------------------------- 
    811   !!   'key_ice_lim'                                     LIM sea-ice model 
    912   !!---------------------------------------------------------------------- 
     13   !!---------------------------------------------------------------------- 
    1014   !!   lim_thd_zdf  : vertical accr./abl. and lateral ablation of sea ice 
    1115   !!---------------------------------------------------------------------- 
    12    !! * Modules used 
    1316   USE par_oce          ! ocean parameters 
    1417   USE phycst           ! ??? 
    … …  
    2225   PRIVATE 
    2326 
    24    !! * Routine accessibility 
    25    PUBLIC lim_thd_zdf        ! called by lim_thd 
    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 
     27   PUBLIC   lim_thd_zdf        ! called by lim_thd 
     28 
     29   REAL(wp) ::   epsi20 = 1.e-20  ,  &  ! constant values 
     30      &          epsi13 = 1.e-13  ,  & 
     31      &          zzero  = 0.e0    ,  & 
     32      &          zone   = 1.e0 
    3333   !!---------------------------------------------------------------------- 
    3434   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    3535   !! $Id$ 
    36    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    37    !!---------------------------------------------------------------------- 
     36   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     37   !!---------------------------------------------------------------------- 
     38 
    3839CONTAINS 
    3940 
    … …  
    6465      !!              - Performs lateral ablation 
    6566      !! 
    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   
     67      !! References : 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   
     69      !!------------------------------------------------------------------ 
     70      INTEGER, INTENT(in) ::   kideb    ! Start point on which the  the computation is applied 
     71      INTEGER, INTENT(in) ::   kiut     ! End point on which the  the computation is applied 
    6972      !! 
    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 
    8073      INTEGER ::   ji       ! dummy loop indices 
    81  
    82       REAL(wp) , DIMENSION(jpij,2) ::  & 
    83          zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 ) 
    84  
     74      REAL(wp), DIMENSION(jpij,2) ::   zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 ) 
    8575      REAL(wp), DIMENSION(jpij) ::  & 
    8676         ztsmlt      &    ! snow/ice surface melting temperature 
    … …  
    9787         , zts_old   &    ! previous surface temperature 
    9888         , zidsn , z1midsn , zidsnic ! tempory variables 
    99  
    100       REAL(wp), DIMENSION(jpij) :: & 
     89      REAL(wp), DIMENSION(jpij) ::   & 
    10190          zfnet       &  ! net heat flux at the top surface( incl. conductive heat flux) 
    10291          , zsprecip  &    ! snow accumulation 
    … …  
    10998          , zfrld_1d    &    ! new sea/ice fraction 
    11099          , zep            ! internal temperature of the 2nd layer of the snow/ice system 
    111  
    112100       REAL(wp), DIMENSION(3) :: &  
    113101          zplediag  &    ! principle diagonal, subdiag. and supdiag. of the  
    … …  
    115103          , zsupdiag  &    ! of the temperatures inside the snow-ice system 
    116104          , zsmbr          ! second member 
    117  
    118105       REAL(wp) :: &  
    119106          zhsu     &     ! thickness of surface layer 
    … …  
    130117          , zb2 , zd2 , zb3 , zd3 & 
    131118          , ztint          ! equivalent temperature at the snow-ice interface 
    132  
    133119       REAL(wp) :: &  
    134120          zexp      &     ! exponential function of the ice thickness 
    … …  
    148134          , zdtic        &  ! ice internal temp. increment 
    149135          , zqnes          ! conductive energy due to ice melting in the first ice layer 
    150  
    151136       REAL(wp) :: &  
    152137          ztbot     &      ! temperature at the bottom surface 
    … …  
    162147          , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables 
    163148          , ztb2, ztb3 
    164  
    165149       REAL(wp) :: &  
    166150          zdrmh         &   ! change in snow/ice thick. after snow-ice formation 
    … …  
    181165       !     Computation of energies due to surface and bottom melting  
    182166       !----------------------------------------------------------------------- 
    183         
    184167        
    185168       DO ji = kideb , kiut 
    … …  
    201184       END DO 
    202185 
    203  
    204186       !------------------------------------------- 
    205187       !  2. Calculate some intermediate variables.   
    … …  
    265247       !     - qstbif_1d, total energy stored in brine pockets (updating) 
    266248       !------------------------------------------------------------------- 
    267  
    268249 
    269250       DO ji = kideb , kiut 
    … …  
    288269       END DO 
    289270 
    290  
    291271       !-------------------------------------------------------------------------------- 
    292272       !  4. Computation of the surface temperature : determined by considering the  
    … …  
    333313          !---computation of the energy balance function  
    334314          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 
     315             &      - qns_ice_1d(ji)                & ! total non solar radiation 
     316             &      - zfcsu (ji)                      ! conductive heat flux from the surface 
    337317          !---computation of surface temperature increment   
    338318          zdts    = -zfts / zdfts 
    … …  
    360340          sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 
    361341#if ! defined key_coupled 
    362           qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
     342          qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
    363343          qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
    364344#endif 
    … …  
    366346       END DO 
    367347 
    368       
    369  
    370348       !     5.2. Calculate available heat for surface ablation.  
    371349       !--------------------------------------------------------------------- 
    372350 
    373351       DO ji = kideb, kiut 
    374           zfnet(ji) = qnsr_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)           
     352          zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)           
    375353          zfnet(ji) = MAX( zzero , zfnet(ji) ) 
    376354          zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) ) 
    … …  
    730708          dvnbq_1d(ji) = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) 
    731709          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 
     710          !---  volume change of ice and snow (used for ocean-ice freshwater flux computation) 
     711          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic 
    735712          rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn 
    736713 
    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  
    742714          !---  Actualize new snow and ice thickness. 
    743  
    744715          h_snow_1d(ji)  = zhsnnew 
    745           h_ice_1d(ji)  = zhicnew 
     716          h_ice_1d (ji)  = zhicnew 
    746717 
    747718       END DO 
    … …  
    799770          qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji) 
    800771          frld_1d(ji)    = zfrld_1d(ji) 
    801  
    802        END DO 
    803         
     772          ! 
     773       END DO 
     774       !  
    804775    END SUBROUTINE lim_thd_zdf 
     776 
    805777#else 
    806    !!====================================================================== 
    807    !!                       ***  MODULE limthd_zdf   *** 
    808    !!                              no sea ice model   
    809    !!====================================================================== 
     778   !!---------------------------------------------------------------------- 
     779   !!   Default Option                                     NO sea-ice model 
     780   !!---------------------------------------------------------------------- 
    810781CONTAINS 
    811782   SUBROUTINE lim_thd_zdf          ! Empty routine 
    812783   END SUBROUTINE lim_thd_zdf 
    813784#endif 
    814  END MODULE limthd_zdf 
     785 
     786   !!====================================================================== 
     787END MODULE limthd_zdf 

  • trunk/NEMO/LIM_SRC/limtrp.F90

    r699 r717  
    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 

  • trunk/NEMO/LIM_SRC/limwri.F90

    r709 r717  
    1515   !!   lim_wri_init : initialization and namelist read 
    1616   !!---------------------------------------------------------------------- 
     17   USE phycst 
    1718   USE dom_oce 
     19   USE daymod 
    1820   USE ice_oce         ! ice variables 
     21   USE sbc_oce 
     22   USE sbc_ice 
    1923   USE dom_ice 
    20    USE sbc_oce         ! surface boundary condition: ocean 
    21    USE phycst 
    22    USE daymod 
    2324   USE ice 
     25 
     26   USE lbclnk 
     27   USE dianam    ! build name of file (routine) 
    2428   USE in_out_manager 
    25    USE dianam    ! build name of file (routine) 
    26    USE lbclnk 
    2729   USE ioipsl 
    2830 
    … …  
    3032   PRIVATE 
    3133 
    32    PUBLIC   lim_wri        ! routine called by lim_step.F90 
     34   PUBLIC   lim_wri        ! routine called by sbcice_lim module 
    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    !! $Id$ 
     53   !! * Substitutions 
     54#   include "vectopt_loop_substitute.h90" 
     55   !!---------------------------------------------------------------------- 
     56   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006) 
     57   !! $Header: $ 
    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  
    … …  
    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 : 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 
     151!!gm        zcmo(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj) 
     152!!gm        zcmo(ji,jj,12) = fsolar (ji,jj) 
     153!!gm        zcmo(ji,jj,13) = fnsolar(ji,jj) 
     154            zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj) 
     155            zcmo(ji,jj,12) = qsr(ji,jj) 
     156            zcmo(ji,jj,13) = qns(ji,jj) 
    157157            ! 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) 
     158!!gm        zcmo(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     159            zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ??? 
     160            zcmo(ji,jj,15) = utaui_ice(ji,jj) 
     161            zcmo(ji,jj,16) = vtaui_ice(ji,jj) 
     162            zcmo(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice(ji,jj) + frld(ji,jj) * qsr(ji,jj)   !!gm a revoir 
     163            zcmo(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qns_ice(ji,jj) + frld(ji,jj) * qns(ji,jj)   !!gm a revoir 
    163164            zcmo(ji,jj,19) = sprecip(ji,jj) 
    164165         END DO 
    … …  
    175176         END DO 
    176177          
    177          IF ( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN  
     178         IF( jf == 7  .OR. jf == 8  .OR. jf == 11 .OR. jf == 12 .OR. jf == 15 .OR.   & 
     179            &                            jf == 23 .OR. jf == 24 .OR. jf == 16 ) THEN  
    178180            CALL lbc_lnk( zfield, 'T', -1. ) 
    179181         ELSE  
    … …  
    181183         ENDIF 
    182184          
    183          IF ( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 
     185         IF( nc(jf) == 1 )  CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 
    184186          
    185187      END DO 
    186188       
    187       IF ( ( nfice * niter + nit000 - 1 ) >= nitend ) THEN 
    188          CALL histclo( nice )  
    189       ENDIF 
     189      IF( ( nn_fsbc * niter + nit000 - 1 ) >= nitend )   CALL histclo( nice )  
    190190      ! 
    191191   END SUBROUTINE lim_wri 
    … …  
    225225         field_13, field_14, field_15, field_16, field_17, field_18,   & 
    226226         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 ) 
     227      !!------------------------------------------------------------------- 
     228 
     229      REWIND ( numnam_ice )                ! Read Namelist namicewri 
    236230      READ   ( numnam_ice  , namiceout ) 
    237231       
    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 
     232      zfield( 1) = field_1 
     233      zfield( 2) = field_2 
     234      zfield( 3) = field_3 
     235      zfield( 4) = field_4 
     236      zfield( 5) = field_5 
     237      zfield( 6) = field_6 
     238      zfield( 7) = field_7 
     239      zfield( 8) = field_8 
     240      zfield( 9) = field_9 
    247241      zfield(10) = field_10 
    248242      zfield(11) = field_11 
    … …  
    274268         DO nf = 1 , noumef          
    275269            WRITE(numout,*) '   ', titn(nf), '   ', nam(nf),'      ', uni(nf),'  ', nc(nf),'        ', cmulti(nf),   & 
    276                '        ', cadd(nf) 
     270               &       '        ', cadd(nf) 
    277271         END DO 
    278272      ENDIF 

  • trunk/NEMO/LIM_SRC/limwri_dimg.h90

    r699 r717  
    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 
    115115          zcmo(ji,jj,9)  = sst_io(ji,jj) 
    … …  
    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 
    158158                rcmoy(ji,jj,9)  = sst_io(ji,jj) 
    … …  
    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 

  • trunk/NEMO/LIM_SRC/thd_ice.F90

    r699 r717  
    5656      fr1_i0_1d   ,     &  !:    "                  "      fr1_i0 
    5757      fr2_i0_1d   ,     &  !:    "                  "      fr2_i0 
    58       qnsr_ice_1d ,     &  !:    "                  "      qns_ice 
     58      qns_ice_1d ,     &  !:    "                  "      qns_ice 
    5959      qfvbq_1d    ,     &  !:    "                  "      qfvbq 
    6060      sist_1d     ,     &  !:    "                  "      sist 

  • trunk/NEMO/OPA_SRC/DIA/diawri.F90

    r714 r717  
    1616   USE sbc_oce         ! surface boundary condition: ocean 
    1717   USE sbc_ice         ! surface boundary condition: ice 
     18   USE sbcssr          ! restoring term toward SST/SSS climatology 
    1819   USE phycst          ! physical constants 
    1920   USE ocfzpt          ! ocean freezing point 
    … …  
    243244            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    244245#endif 
    245 #if ! defined key_dynspg_rl && defined key_ice_lim 
    246          ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to 
    247          !    internal damping to Levitus that can be diagnosed from others 
    248          ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup 
    249          CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater"          , "kg/m2/s",   &  ! fsalt 
    250             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    251          CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater"        , "kg/m2/s",   &  ! fmass 
    252             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    253 #endif 
     246!!$#if ! defined key_dynspg_rl && defined key_ice_lim 
     247!!$         ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to 
     248!!$         !    internal damping to Levitus that can be diagnosed from others 
     249!!$         ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup 
     250!!$         CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater"          , "kg/m2/s",   &  ! fsalt 
     251!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     252!!$         CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater"        , "kg/m2/s",   &  ! fmass 
     253!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     254!!$#endif 
    254255         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! emp 
    255256            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    256          CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs 
    257             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     257!!$         CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs 
     258!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    258259         CALL histdef( nid_T, "sowaflcd", "concentration/dilution water flux"  , "kg/m2/s",   &  ! emps 
    259260            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    … …  
    421422      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height 
    422423#endif 
    423 #if ! defined key_dynspg_rl && defined key_ice_lim 
    424       CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:)    , ndim_hT, ndex_hT )   ! ice=>ocean water flux 
    425       CALL histwrite( nid_T, "sowaflep", it, fmass(:,:)    , ndim_hT, ndex_hT )   ! atmos=>ocean water flux 
    426 #endif 
     424!!$#if ! defined key_dynspg_rl && defined key_ice_lim 
     425!!$      CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:)    , ndim_hT, ndex_hT )   ! ice=>ocean water flux 
     426!!$      CALL histwrite( nid_T, "sowaflep", it, fmass(:,:)    , ndim_hT, ndex_hT )   ! atmos=>ocean water flux 
     427!!$#endif 
    427428      CALL histwrite( nid_T, "sowaflup", it, emp           , ndim_hT, ndex_hT )   ! upward water flux 
    428       CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff 
     429!!$      CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff 
    429430      CALL histwrite( nid_T, "sowaflcd", it, emps          , ndim_hT, ndex_hT )   ! c/d water flux 
    430431      zw2d(:,:) = emps(:,:) * sn(:,:,1) * tmask(:,:,1) 

  • trunk/NEMO/OPA_SRC/DIA/diawri_dimg.h90

    r714 r717  
    7676    !! * modules used 
    7777    USE lib_mpp 
    78     USE dtasst, ONLY : sst 
    7978 
    8079    !! * Arguments 

  • trunk/NEMO/OPA_SRC/DOM/domain.F90

    r709 r717  
    143143      NAMELIST/namrun/ no    , cexper   , ln_rstart , nrstdt , nit000,         & 
    144144         &             nitend, ndate0   , nleapy   , ninist , nstock,          & 
    145          &             nwrite, nrunoff  , ln_dimgnnn 
     145         &             nwrite, ln_dimgnnn 
    146146 
    147147      NAMELIST/namdom/ ntopo , e3zps_min, e3zps_rat, ngrid  , nmsh  ,   & 
    148148         &             nacc  , atfp     , rdt      , rdtmin , rdtmax,   & 
    149          &             rdth  , rdtbt    , nfice    , nfbulk , nclosea 
     149         &             rdth  , rdtbt    , nclosea 
    150150      NAMELIST/namcla/ n_cla 
    151151      !!---------------------------------------------------------------------- 
    … …  
    174174         WRITE(numout,*) '           frequency of restart file       nstock    = ', nstock 
    175175         WRITE(numout,*) '           frequency of output file        nwrite    = ', nwrite 
    176          WRITE(numout,*) '           runoff option                   nrunoff   = ', nrunoff 
    177176         WRITE(numout,*) '           multi file dimgout           ln_dimgnnn   = ', ln_dimgnnn 
    178177      ENDIF 
    … …  
    258257      ENDIF 
    259258 
    260       IF( lk_ice_lim ) THEN 
    261          IF(lwp) WRITE(numout,*) '           ice model coupling frequency      nfice  = ', nfice 
    262          nfbulk = nfice 
    263          IF( MOD( rday, nfice*rdt ) /= 0 ) THEN  
    264             IF(lwp) WRITE(numout,*) ' ' 
    265             IF(lwp) WRITE(numout,*) 'W A R N I N G :  nfice is NOT a multiple of the number of time steps in a day' 
    266             IF(lwp) WRITE(numout,*) ' ' 
    267          ENDIF 
    268          IF(lwp) WRITE(numout,*) '           bulk computation frequency       nfbulk  = ', nfbulk, ' = nfice if ice model used' 
    269          IF(lwp) WRITE(numout,*) '           flag closed sea or not           nclosea = ', nclosea 
    270       ENDIF 
    271  
    272259      ! Default values 
    273260      n_cla = 0 

  • trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90

    r708 r717  
    1717   USE dom_oce         ! ocean space and time domain 
    1818   USE sbc_oce         ! surface boundary condition: ocean 
     19   USE sbcrnf          ! river runoffs 
    1920   USE ldftra_oce      ! ocean tracer   lateral physics 
    2021   USE phycst          ! physical constants 
    … …  
    177178            DO ji = 1, jpi 
    178179               zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  & 
    179                   &      + aht0 * upsrnfh(ji,jj)                          ! enhanced near river mouths 
     180                  &      + aht0 * rnfmsk(ji,jj)                          ! enhanced near river mouths 
    180181               ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 ) 
    181182               ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 ) 
    … …  
    352353            DO ji = 1, jpi 
    353354               zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  & 
    354                   &      + aht0 * upsrnfh(ji,jj)                          ! enhanced near river mouths 
     355                  &      + aht0 * rnfmsk(ji,jj)                          ! enhanced near river mouths 
    355356               ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 ) 
    356357               ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 ) 

  • trunk/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r713 r717  
    1818   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns       !: sea heat flux: non solar                     [W/m2] 
    1919   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr       !: sea heat flux:     solar                     [W/m2] 
    20    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qrp       !: heat flux damping                            [w/m2] 
    2120   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp       !: freshwater budget: volume flux               [Kg/m2/s] 
    2221   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps      !: freshwater budget: concentration/dillution   [Kg/m2/s] 
    23    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   erp       !: evaporation damping                          [kg/m2/s] 
    2422 
    2523   !!---------------------------------------------------------------------- 

  • trunk/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r702 r717  
    104104         DO jj = 2, jpj 
    105105            DO ji = fs_2, jpi   ! vector opt. 
    106                u_oce(ji,jj) = 0.5 * ( ssu_m(ji-1,jj  ) + ssu_m(ji-1,jj-1) ) 
    107                v_oce(ji,jj) = 0.5 * ( ssv_m(ji  ,jj-1) + ssv_m(ji-1,jj-1) ) 
     106               ui_oce(ji,jj) = 0.5 * ( ssu_m(ji-1,jj  ) + ssu_m(ji-1,jj-1) ) 
     107               vi_oce(ji,jj) = 0.5 * ( ssv_m(ji  ,jj-1) + ssv_m(ji-1,jj-1) ) 
    108108            END DO 
    109109         END DO 
    110          CALL lbc_lnk( u_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices) 
    111          CALL lbc_lnk( v_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices) 
     110         CALL lbc_lnk( ui_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices) 
     111         CALL lbc_lnk( vi_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices) 
    112112 
    113113         ! ... masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
    … …  
    124124            CALL blk_ice_clio() 
    125125         CASE( 4 )           ! CORE bulk formulation 
    126             CALL blk_ice_core( sist     , u_ice    , v_ice    , alb_ice_cs,                      & 
     126            CALL blk_ice_core( sist     , ui_ice   , vi_ice   , alb_ice_cs,                      & 
    127127               &                                     utaui_ice, vtaui_ice , qns_ice , qsr_ice,   & 
    128128               &                                     qla_ice  , dqns_ice  , dqla_ice,            & 
    … …  
    148148         ;                              CALL lim_thd      ( kt )      ! Ice thermodynamics  
    149149         ;                              CALL lim_sbc      ( kt )      ! Ice/Ocean Mass & Heat fluxes  
    150          IF( MOD( kt+nfice-1, ninfo ) == 0 .OR.   & 
     150         IF( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR.   & 
    151151            &  ntmoy == 1 )             CALL lim_dia      ( kt )      ! Ice Diagnostics  
    152152         ;                              CALL lim_wri      ( kt )      ! Ice outputs  

  • trunk/NEMO/OPA_SRC/SBC/sbcssr.F90

    r713 r717  
    3333   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sst   ! structure of input SST (file informations, fields read) 
    3434   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sss   ! structure of input SSS (file informations, fields read) 
     35 
     36   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   erp       !: evaporation damping                          [kg/m2/s] 
     37   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qrp       !: heat flux damping                            [w/m2] 
    3538 
    3639   !! * Namelist namsbc_ssr 
    … …  
    99102            WRITE(numout,*) '                           (Yes=2, volume flux) ' 
    100103            WRITE(numout,*) '          dQ/dT (restoring magnitude on SST)     dqdt     = ', dqdt, ' W/m2/K' 
    101             WRITE(numout,*) '          dE/dS (restoring magnitude on SST)     deds     = ', deds, ' ???' 
     104            WRITE(numout,*) '          dE/dS (restoring magnitude on SST)     deds     = ', deds, ' ...' 
    102105         ENDIF 
    103106 

  • trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90

    r708 r717  
    11MODULE traadv_cen2 
    2    !!============================================================================== 
    3    !!                       ***  MODULE  traadv_cen2  *** 
     2   !!====================================================================== 
     3   !!                     ***  MODULE  traadv_cen2  *** 
    44   !! Ocean active tracers:  horizontal & vertical advective trend 
    5    !!============================================================================== 
    6    !! History :  8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv  
    7    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
    8    !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
    9    !!            " "  !  06-04  (R. Benshila, G. Madec) Step reorganization 
     5   !!====================================================================== 
     6   !! History :   8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad=traadv  
     7   !!             8.5  !  02-06  (G. Madec)  F90: Free form and module 
     8   !!             9.0  !  04-08  (C. Talandier) New trends organization 
     9   !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization 
     10   !!             " "  !  06-04  (R. Benshila, G. Madec) Step reorganization 
     11   !!             " "  !  06-07  (G. madec)  add ups_orca_set routine 
    1012   !!---------------------------------------------------------------------- 
    1113 
    … …  
    1416   !!                  vertical advection trends using a seconder order 
    1517   !!                  centered scheme. (k-j-i loops) 
     18   !!   ups_orca_set : allow mixed upstream/centered scheme in specific 
     19   !!                  area (set for orca 2 and 4 only) 
    1620   !!---------------------------------------------------------------------- 
    1721   USE oce             ! ocean dynamics and active tracers 
    … …  
    2125   USE trdmod_oce      ! ocean variables trends 
    2226   USE trdmod          ! ocean active tracers trends  
     27   USE closea          ! closed sea 
    2328   USE trabbl          ! advective term in the BBL 
    2429   USE ocfzpt          ! 
     30   USE sbcrnf          ! river runoffs 
    2531   USE in_out_manager  ! I/O manager 
    2632   USE lib_mpp 
    … …  
    3238   PRIVATE 
    3339 
    34    PUBLIC   tra_adv_cen2   ! routine called by step.F90 
     40   PUBLIC   tra_adv_cen2    ! routine called by step.F90 
     41   PUBLIC   ups_orca_set    ! routine used by traadv_cen2_jki.F90 
     42 
     43   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   upsmsk    !: mixed upstream/centered scheme near some straits  
     44   !                                                   !  and in closed seas (orca 2 and 4 configurations) 
    3545 
    3646   REAL(wp), DIMENSION(jpi,jpj) ::   btr2   ! inverse of T-point surface [1/(e1t*e2t)] 
    … …  
    4050#  include "vectopt_loop_substitute.h90" 
    4151   !!---------------------------------------------------------------------- 
    42    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     52   !!   OPA 9.0 , LOCEAN-IPSL (2006)  
    4353   !! $Id$ 
    4454   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    119129      !! 
    120130      INTEGER  ::   ji, jj, jk                           ! dummy loop indices 
    121       REAL(wp) ::                           & 
    122          zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars 
    123          zhw, ze3tr, zcofi, zcofj,          &  !    "         " 
    124          zupsut, zupsvt, zupsus, zupsvs,    &  !    "         " 
    125          zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         " 
    126          zcofk, zupst, zupss, zcent,        &  !    "         " 
    127          zcens, zfp_w, zfm_w,               &  !    "         " 
    128          zcenut, zcenvt, zcenus, zcenvs,    &  !    "         " 
    129          z_hdivn_x, z_hdivn_y, z_hdivn 
    130       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind   ! 3D workspace  
     131      REAL(wp) ::   zta, zsa, zbtr, zhw, ze3tr,       &  ! temporary scalars 
     132         &          zfp_ui, zfp_vj, zfp_w , zfui  ,   &  !    "         " 
     133         &          zfm_ui, zfm_vj, zfm_w , zfvj  ,   &  !    "         " 
     134         &          zcofi , zcofj , zcofk ,           &  !    "         " 
     135         &          zupsut, zupsus, zcenut, zcenus,   &  !    "         " 
     136         &          zupsvt, zupsvs, zcenvt, zcenvs,   &  !    "         " 
     137         &          zupst , zupss , zcent , zcens ,   &  !    "         " 
     138         &          z_hdivn_x, z_hdivn_y, z_hdivn  
     139      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind   ! 3D workspace 
    131140      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww, ztrds         !  "      " 
    132141      !!---------------------------------------------------------------------- 
    … …  
    137146         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case' 
    138147         IF(lwp) WRITE(numout,*) 
    139          !  
    140          btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
     148         ! 
     149         upsmsk(:,:) = 0.e0                              ! not upstream by default 
     150         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits 
     151         !                                               ! and in closed seas (orca2 and orca4 only) 
     152         !    
     153         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )        ! inverse of T-point surface 
    141154      ENDIF 
    142155 
    … …  
    146159         DO jj = 1, jpj 
    147160            DO ji = 1, jpi 
    148                zind(ji,jj,jk) =  MAX ( upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff 
    149                   &                    upsadv(ji,jj)                     &  ! in the vicinity of some straits 
     161               zind(ji,jj,jk) = MAX (   & 
     162                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows) 
     163                  upsmsk(ji,jj)                      &  ! some of some straits 
    150164#if defined key_ice_lim 
    151                   &                  , tmask(ji,jj,jk)                   &  ! half upstream tracer fluxes 
    152                   &                  * MAX( 0., SIGN( 1., fzptn(ji,jj)   &  ! if tn < ("freezing"+0.1 ) 
    153                   &                                +0.1-tn(ji,jj,jk) ) ) & 
     165                  !                                     ! below ice covered area (if tn < "freezing"+0.1 ) 
     166                  , MAX(  0., SIGN( 1., fzptn(ji,jj) + 0.1 - tn(ji,jj,jk) )  ) * tmask(ji,jj,jk)   & 
    154167#endif 
    155168                  &                  ) 
    … …  
    158171      END DO 
    159172 
    160  
    161       !  Horizontal advective fluxes 
    162       ! ----------------------------- 
     173      ! I. Horizontal advective fluxes 
     174      ! ------------------------------ 
     175      !  Second order centered tracer flux at u and v-points 
     176      ! ----------------------------------------------------- 
    163177      !                                                ! =============== 
    164178      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    208222               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 
    209223#endif 
    210                ! horizontal advective trends 
    211                zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
     224               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &    ! horizontal advective trends 
    212225                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    213226               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   & 
    214227                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  ) 
    215                ! add it to the general tracer trends 
    216                ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
     228               ta(ji,jj,jk) = ta(ji,jj,jk) + zta                          ! add it to the general tracer trends 
    217229               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
    218230            END DO 
    … …  
    279291         &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    280292 
    281       ! 4. "zonal" mean advective heat and salt transport  
    282       ! ------------------------------------------------- 
    283  
     293      ! "zonal" mean advective heat and salt transport 
    284294      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
    285295         IF( lk_zco ) THEN 
    … …  
    313323      ENDIF 
    314324 
    315       ! 1. Vertical advective fluxes 
     325      ! 1. Vertical advective fluxes (Second order centered tracer flux at w-point) 
    316326      ! ---------------------------- 
    317       ! Second order centered tracer flux at w-point 
    318327      DO jk = 2, jpk 
    319328         DO jj = 2, jpjm1 
    320329            DO ji = fs_2, fs_jpim1   ! vector opt. 
    321                ! upstream indicator 
    322                zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
    323                ! velocity * 1/2 
    324                zhw = 0.5 * pwn(ji,jj,jk) 
    325                ! upstream scheme 
    326                zfp_w = zhw + ABS( zhw ) 
     330               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )         ! upstream indicator 
     331               zhw = 0.5 * pwn(ji,jj,jk)                               ! velocity * 1/2 
     332               zfp_w = zhw + ABS( zhw )                                ! upstream scheme 
    327333               zfm_w = zhw - ABS( zhw ) 
    328334               zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1) 
    329335               zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1) 
    330                ! centered scheme 
    331                zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) 
     336               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )         ! centered scheme 
    332337               zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) 
    333                ! mixed centered / upstream scheme 
    334                zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent 
     338               zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent      ! mixed centered / upstream scheme 
    335339               zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens 
    336340            END DO 
    … …  
    344348            DO ji = fs_2, fs_jpim1   ! vector opt. 
    345349               ze3tr = 1. / fse3t(ji,jj,jk) 
    346                ! vertical advective trends 
    347                zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
     350               zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )     ! vertical advective trends 
    348351               zsa = - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) 
    349                ! add it to the general tracer trends 
    350                ta(ji,jj,jk) =  ta(ji,jj,jk) + zta 
     352               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta                      ! add it to the general tracer trends 
    351353               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa 
    352354            END DO 
    … …  
    387389      ! 
    388390   END SUBROUTINE tra_adv_cen2 
     391    
     392    
     393   SUBROUTINE ups_orca_set 
     394      !!---------------------------------------------------------------------- 
     395      !!                  ***  ROUTINE ups_orca_set  *** 
     396      !!        
     397      !! ** Purpose :   add a portion of upstream scheme in area where the 
     398      !!                centered scheme generates too strong overshoot 
     399      !! 
     400      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk 
     401      !!                array to nozero value in some straith.  
     402      !! 
     403      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca 
     404      !!---------------------------------------------------------------------- 
     405      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers 
     406      !!---------------------------------------------------------------------- 
     407       
     408      ! mixed upstream/centered scheme near river mouths 
     409      ! ------------------------------------------------ 
     410      SELECT CASE ( jp_cfg ) 
     411      !                                        ! ======================= 
     412      CASE ( 4 )                               !  ORCA_R4 configuration  
     413         !                                     ! ======================= 
     414         !                                          ! Gibraltar Strait 
     415         ii0 =  70   ;   ii1 =  71 
     416         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 
     417         ! 
     418         !                                     ! ======================= 
     419      CASE ( 2 )                               !  ORCA_R2 configuration  
     420         !                                     ! ======================= 
     421         !                                          ! Gibraltar Strait 
     422         ij0 = 102   ;   ij1 = 102 
     423         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20 
     424         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40 
     425         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 
     426         ij0 = 101   ;   ij1 = 102 
     427         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 
     428         !                                          ! Bab el Mandeb Strait 
     429         ij0 =  87   ;   ij1 =  88 
     430         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10 
     431         ij0 =  88   ;   ij1 =  88 
     432         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 
     433         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40 
     434         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 
     435         ij0 =  89   ;   ij1 =  89 
     436         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 
     437         ij0 =  90   ;   ij1 =  90 
     438         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 
     439         !                                          ! Sound Strait 
     440         ij0 = 116   ;   ij1 = 116 
     441         ii0 = 144   ;   ii1 = 144   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 
     442         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 
     443         ii0 = 148   ;   ii1 = 148   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 
     444         ! 
     445      END SELECT  
     446       
     447      ! mixed upstream/centered scheme over closed seas 
     448      ! ----------------------------------------------- 
     449      CALL clo_ups( upsmsk(:,:) ) 
     450      ! 
     451   END SUBROUTINE ups_orca_set 
    389452 
    390453   !!====================================================================== 

  • trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90

    r708 r717  
    44   !! Ocean active tracers:  horizontal & vertical advective trend 
    55   !!============================================================================== 
    6    !! History :  8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv 
    7    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
    8    !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
    9    !!            " "  !  06-04  (R. Benshila, G. Madec) Step reorganization 
     6   !! History :   8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv 
     7   !!             8.5  !  02-06  (G. Madec)  F90: Free form and module 
     8   !!             9.0  !  04-08  (C. Talandier) New trends organization 
     9   !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization 
     10   !!             " "  !  06-04  (R. Benshila, G. Madec) Step reorganization 
     11   !!             " "  !  06-07  (G. madec)  add ups_orca_set routine 
    1012   !!---------------------------------------------------------------------- 
    1113 
    1214   !!---------------------------------------------------------------------- 
    1315   !!   tra_adv_cen2_jki : update the tracer trend with the horizontal and 
    14    !!                  vertical advection trends using a seconder order 
    15    !!                  centered scheme. Auto-tasking case, k-slab for 
    16    !!                  hor. adv., j-slab for vert. adv. (j-k-i loops) 
     16   !!                      vertical advection trends using a seconder order 
     17   !!                      centered scheme. Auto-tasking case, k-slab for 
     18   !!                      hor. adv., j-slab for vert. adv. (j-k-i loops) 
    1719   !!---------------------------------------------------------------------- 
    1820   USE oce             ! ocean dynamics and active tracers 
    … …  
    2628   USE lib_mpp 
    2729   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
     30   USE sbcrnf          ! river runoffs 
    2831   USE in_out_manager  ! I/O manager 
    2932   USE diaptr          ! poleward transport diagnostics 
    3033   USE prtctl          ! Print control 
    31  
     34    
    3235   IMPLICIT NONE 
    3336   PRIVATE 
    3437 
    35    PUBLIC   tra_adv_cen2_jki   ! routine called by step.F90 
    36  
    37    REAL(wp), DIMENSION(jpi,jpj) ::   btr2   ! inverse of T-point surface [1/(e1t*e2t)] 
     38   PUBLIC   tra_adv_cen2_jki    ! routine called by step.F90 
     39 
     40   REAL(wp), DIMENSION(jpi,jpj) ::   btr2    ! inverse of T-point surface [1/(e1t*e2t)] 
    3841 
    3942   !! * Substitutions 
    … …  
    8790      !!         zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[tn] ) 
    8891      !!                                      +mj-1( e1v*e3v  vn  mj[tn] )  } 
     92      !!         C A U T I O N : the trend saved is the centered trend only. 
     93      !!      It does not take into account the upstream part of the scheme. 
    8994      !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so  
    9095      !!      they vanish from the expression of the flux and divergence. 
    … …  
    109114      !! 
    110115      !! ** Action :  - update (ta,sa) with the now advective tracer trends 
    111       !!              - save trends in (ztrdt,ztrds) ('key_trdtra') 
    112116      !!---------------------------------------------------------------------- 
    113       USE oce, ONLY :   zwx => ua   ! use ua as workspace 
    114       USE oce, ONLY :   zwy => va   ! use va as workspace 
    115       !! 
    116       INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index 
    117       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component 
    118       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component 
    119       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component 
    120       !! 
     117      USE oce               ,   zwx => ua          ! use ua as workspace 
     118      USE oce               ,   zwy => va          ! use va as workspace 
     119      USE traadv_cen2, ONLY :   ups_orca_set,   &  ! upstream indicator near some straits  
     120         &                      upsmsk             ! and over closed sea (orca 2 and 4) 
     121 
     122      INTEGER , INTENT(in)                         ::   kt              ! ocean time-step index 
     123      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! now ocean velocity fields 
     124 
    121125      INTEGER  ::   ji, jj, jk                           ! dummy loop indices 
    122       REAL(wp) ::                           & 
    123          zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars 
    124          zhw, ze3tr, zcofi, zcofj,          &  !    "         " 
    125          zupsut, zupsvt, zupsus, zupsvs,    &  !    "         " 
    126          zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         " 
    127          zcofk, zupst, zupss, zcent,        &  !    "         " 
    128          zcens, zfp_w, zfm_w,               &  !    "         " 
    129          zcenut, zcenvt, zcenus, zcenvs,    &  !    "         " 
    130          z_hdivn_x, z_hdivn_y, z_hdivn 
    131       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind   ! 3D workspace  
    132       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww, ztrds         !  "      " 
     126      REAL(wp) ::   zta, zsa, zbtr, zhw, ze3tr,       &  ! temporary scalars 
     127         &          zfp_ui, zfp_vj, zfp_w , zfui  ,   &  !    "         " 
     128         &          zfm_ui, zfm_vj, zfm_w , zfvj  ,   &  !    "         " 
     129         &          zcofi , zcofj , zcofk ,           &  !    "         " 
     130         &          zupsut, zupsus, zcenut, zcenus,   &  !    "         " 
     131         &          zupsvt, zupsvs, zcenvt, zcenvs,   &  !    "         " 
     132         &          zupst , zupss , zcent , zcens ,   &  !    "         " 
     133         &          z_hdivn_x, z_hdivn_y, z_hdivn  
     134      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind,   &  ! temporary workspace arrays 
     135         &                                  zww, ztrds             !    "         " 
    133136      !!---------------------------------------------------------------------- 
    134  
     137      ! 
    135138      IF( kt == nit000 ) THEN 
    136139         IF(lwp) WRITE(numout,*) 
    … …  
    138141         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   Auto-tasking case' 
    139142         IF(lwp) WRITE(numout,*) 
    140          ! 
    141          btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
     143         !          
     144         upsmsk(:,:) = 0.e0                              ! not upstream by default 
     145         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits 
     146         !                                               ! and in closed seas (orca2 and orca4 only) 
     147 
     148         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )        ! inverse of T-point surface 
    142149      ENDIF 
    143150 
    … …  
    147154      DO jk = 1, jpkm1                                 ! Horizontal slab 
    148155         !                                             ! =============== 
    149          !  Upstream / centered scheme indicator 
    150          ! -------------------------------------- 
     156 
     157         ! 0. Upstream / centered scheme indicator 
     158         ! --------------------------------------- 
    151159         DO jj = 1, jpj 
    152160            DO ji = 1, jpi 
    153                zind(ji,jj,jk) =  MAX (              & 
    154                   upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff 
    155                   upsadv(ji,jj)                     &  ! in the vicinity of some straits 
     161               zind(ji,jj,jk) = MAX (   & 
     162                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows) 
     163                  upsmsk(ji,jj)                      &  ! some of some straits 
    156164#if defined key_ice_lim 
    157                   , tmask(ji,jj,jk)                 &  ! half upstream tracer fluxes if tn < ("freezing"+0.1 ) 
    158                   * MAX( 0., SIGN( 1., fzptn(ji,jj)+0.1-tn(ji,jj,jk) ) )   & 
    159 #endif 
    160                ) 
    161             END DO 
    162          END DO 
    163  
    164          !  Horizontal advective fluxes 
    165          ! ----------------------------- 
     165                  !                                     ! below ice covered area (if tn < "freezing"+0.1 ) 
     166                  , MAX(  0., SIGN( 1., fzptn(ji,jj) + 0.1 - tn(ji,jj,jk) )  ) * tmask(ji,jj,jk)   & 
     167#endif 
     168                  &                  ) 
     169            END DO 
     170         END DO 
     171 
     172 
     173         ! I. Horizontal advective fluxes 
     174         ! ------------------------------ 
    166175         ! Second order centered tracer flux at u and v-points 
    167176         DO jj = 1, jpjm1 
    168177            DO ji = 1, fs_jpim1   ! vector opt. 
    169                ! upstream indicator 
    170                zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) ) 
     178               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )               ! upstream indicator 
    171179               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) ) 
    172                ! volume fluxes * 1/2 
    173 #if defined key_zco 
    174                zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk) 
     180#if defined key_zco 
     181               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)                       ! volume fluxes * 1/2 (zco) 
    175182               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk) 
    176183#else 
    177                zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
     184               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)     ! volume fluxes * 1/2 
    178185               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    179186#endif 
    180                ! upstream scheme 
    181                zfp_ui = zfui + ABS( zfui ) 
     187               zfp_ui = zfui + ABS( zfui )                                   ! upstream scheme 
    182188               zfp_vj = zfvj + ABS( zfvj ) 
    183189               zfm_ui = zfui - ABS( zfui ) 
    … …  
    187193               zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk) 
    188194               zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk) 
    189                ! centered scheme 
    190                zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) ) 
     195               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )           ! centered scheme 
    191196               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) ) 
    192197               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) ) 
    193198               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) ) 
    194                ! mixed centered / upstream scheme 
    195                zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut 
     199               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut          ! mixed centered / upstream scheme 
    196200               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt 
    197201               zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus 
    … …  
    200204         END DO 
    201205 
    202          ! Tracer flux divergence at t-point added to the general trend 
     206         ! 2. Tracer flux divergence at t-point added to the general trend 
     207         ! --------------------------------------------------------------- 
     208 
    203209         DO jj = 2, jpjm1 
    204210            DO ji = fs_2, fs_jpim1   ! vector opt. 
    205211#if defined key_zco 
    206                zbtr = btr2(ji,jj)  
    207 #else 
    208                zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 
    209 #endif 
    210                ! horizontal advective trends 
    211                zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
     212               zbtr = btr2(ji,jj)                                           ! inverse of the volume (zco) 
     213#else 
     214               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)                         ! inverse of the volume 
     215#endif 
     216               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &      ! horizontal advective trends 
    212217                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    213218               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   & 
    214219                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  ) 
    215                ! add it to the general tracer trends 
    216                ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
     220               ta(ji,jj,jk) = ta(ji,jj,jk) + zta                            ! add it to the general tracer trends 
    217221               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
    218222            END DO 
    … …  
    225229         &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    226230 
    227       ! Save the horizontal advective trends for diagnostics 
     231      ! 3. Save the horizontal advective trends for diagnostic 
    228232      IF( l_trdtra )   THEN 
    229233         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0 
    … …  
    298302      ENDIF 
    299303 
    300       ! Vertical advection 
    301       ! -------------------- 
     304      ! II. Vertical advection 
     305      ! ---------------------- 
    302306!CDIR PARALLEL DO 
    303307!$OMP PARALLEL DO 
    … …  
    312316         IF( lk_dynspg_rl ) THEN                          ! rigid lid : flux set to zero 
    313317            zwz(:,jj, 1 ) = 0.e0    ;    zww(:,jj, 1 ) = 0.e0 
    314          ELSE                                             ! free surface-constant volume 
     318         ELSE                                             ! free surface-constant volume : advection across the surface 
    315319            zwz(:,jj, 1 ) = pwn(:,jj,1) * tn(:,jj,1) 
    316320            zww(:,jj, 1 ) = pwn(:,jj,1) * sn(:,jj,1) 
    317321         ENDIF 
    318322          
    319          ! Vertical advective fluxes at w-point 
     323         ! 1. Vertical advective fluxes (Second order centered tracer flux at w-point) 
     324         ! ---------------------------- 
    320325         DO jk = 2, jpk 
    321326            DO ji = 2, jpim1 
    322                ! upstream indicator 
    323                zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
    324                ! velocity * 1/2 
    325                zhw = 0.5 * pwn(ji,jj,jk) 
    326                ! upstream scheme 
    327                zfp_w = zhw + ABS( zhw ) 
     327               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )          ! upstream indicator 
     328               zhw = 0.5 * pwn(ji,jj,jk)                                ! velocity * 1/2 
     329               zfp_w = zhw + ABS( zhw )                                 ! upstream scheme 
    328330               zfm_w = zhw - ABS( zhw ) 
    329331               zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1) 
    330332               zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1) 
    331                ! centered scheme 
    332                zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) 
     333               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )          ! centered scheme 
    333334               zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) 
    334                ! mixed centered / upstream scheme 
    335                zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent 
     335               zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent       ! mixed centered / upstream scheme 
    336336               zww(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens 
    337337            END DO 
    338338         END DO 
    339339 
    340          ! Tracer flux divergence at t-point added to the general trend 
     340         ! 2. Tracer flux divergence at t-point added to the general trend 
     341         ! ------------------------- 
    341342         DO jk = 1, jpkm1 
    342343            DO ji = 2, jpim1 
    343344               ze3tr = 1. / fse3t(ji,jj,jk) 
    344                ! vertical advective trends 
    345                zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
     345               zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )      ! vertical advective trends 
    346346               zsa = - ze3tr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) 
    347                ! add it to the general tracer trends 
    348                ta(ji,jj,jk) =  ta(ji,jj,jk) + zta 
     347               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta                       ! add it to the general tracer trends 
    349348               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa 
    350349            END DO 
    … …  
    354353      !                                                ! =============== 
    355354 
    356       ! Save the vertical advective trends for diagnostic 
     355      ! 3. Save the vertical advective trends for diagnostic 
     356      ! ---------------------------------------------------- 
    357357      IF( l_trdtra )   THEN 
    358358         ! Recompute the vertical advection zta & zsa trends computed  

  • trunk/NEMO/OPA_SRC/ice_oce.F90

    r709 r717  
    3939      fcalving            !: Iceberg calving  
    4040# endif 
    41  
    42    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: field exchanges with ice model to ocean 
    43       sst_io, sss_io , &  !: sea surface temperature (C) and salinity (PSU) 
    44       u_io  , v_io   , &  !: velocity at ice surface (m/s) 
    45       fsolar, fnsolar, &  !: solar and non-solar heat fluxes (W/m2) 
    46       fsalt , fmass  , &  !: salt and freshwater fluxes 
    47       ftaux , ftauy  , &  !: wind stresses 
    48       gtaux , gtauy       !: wind stresses 
    4941    
    5042   REAL(wp), PUBLIC ::   &  !: 
    … …  
    5951#endif 
    6052 
    61    INTEGER, PUBLIC ::   &  !: namdom : space/time domain (namlist) 
    62       nfice =  5           !: coupling frequency OPA ICELLN  nfice  
    63  
    6453   !!---------------------------------------------------------------------- 
    6554END MODULE ice_oce 

  • trunk/NEMO/OPA_SRC/lbclnk.F90

    r699 r717  
    329329 
    330330 
    331    SUBROUTINE lbc_lnk_3d( pt3d, cd_type, psgn, cd_mpp ) 
     331   SUBROUTINE lbc_lnk_3d( pt3d, cd_type, psgn, cd_mpp, pval ) 
    332332      !!--------------------------------------------------------------------- 
    333333      !!                  ***  ROUTINE lbc_lnk_3d  *** 
    … …  
    355355      CHARACTER(len=3), INTENT( in ), OPTIONAL ::    & 
    356356         cd_mpp        ! fill the overlap area only (here do nothing) 
     357      REAL(wp)        , INTENT(in   ), OPTIONAL           ::   pval      ! background value (used at closed boundaries) 
    357358 
    358359      !! * Local declarations 
    359360      INTEGER  ::   ji, jk 
    360361      INTEGER  ::   ijt, iju 
     362      REAL(wp) ::   zland 
    361363      !!---------------------------------------------------------------------- 
    362364      !!  OPA 9.0 , LOCEAN-IPSL (2005)  
    … …  
    365367      !!---------------------------------------------------------------------- 
    366368 
    367       IF (PRESENT(cd_mpp)) THEN 
     369      IF( PRESENT( pval ) ) THEN      ! set land value (zero by default) 
     370         zland = pval 
     371      ELSE 
     372         zland = 0.e0 
     373      ENDIF 
     374 
     375 
     376      IF( PRESENT( cd_mpp ) ) THEN 
    368377         ! only fill the overlap area and extra allows  
    369378         ! this is in mpp case. In this module, just do nothing 
    … …  
    385394            SELECT CASE ( cd_type ) 
    386395            CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points 
    387                pt3d( 1 ,:,jk) = 0.e0 
    388                pt3d(jpi,:,jk) = 0.e0 
    389             CASE ( 'F' )                               ! F-point 
    390                pt3d(jpi,:,jk) = 0.e0 
     396               pt3d( 1 ,:,jk) = zland 
     397               pt3d(jpi,:,jk) = zland 
     398            CASE ( 'F' )                               ! F-point 
     399               pt3d(jpi,:,jk) = zland 
    391400            END SELECT 
    392401 
    … …  
    402411            CASE ( 'T' , 'U' , 'W' )                   ! T-, U-, W-points 
    403412               pt3d(:, 1 ,jk) = pt3d(:,3,jk) 
    404                pt3d(:,jpj,jk) = 0.e0 
     413               pt3d(:,jpj,jk) = zland 
    405414            CASE ( 'V' , 'F' )                         ! V-, F-points 
    406415               pt3d(:, 1 ,jk) = psgn * pt3d(:,2,jk) 
    407                pt3d(:,jpj,jk) = 0.e0 
     416               pt3d(:,jpj,jk) = zland 
    408417            END SELECT 
    409418 
    410419         CASE ( 3 , 4 )                        ! *  North fold  T-point pivot 
    411420 
    412             pt3d( 1 ,jpj,jk) = 0.e0 
    413             pt3d(jpi,jpj,jk) = 0.e0 
     421            pt3d( 1 ,jpj,jk) = zland 
     422            pt3d(jpi,jpj,jk) = zland 
    414423 
    415424            SELECT CASE ( cd_type ) 
    … …  
    417426               DO ji = 2, jpi 
    418427                  ijt = jpi-ji+2 
    419                   pt3d(ji, 1 ,jk) = 0.e0 
     428                  pt3d(ji, 1 ,jk) = zland 
    420429                  pt3d(ji,jpj,jk) = psgn * pt3d(ijt,jpj-2,jk) 
    421430               END DO 
    … …  
    427436               DO ji = 1, jpi-1 
    428437                  iju = jpi-ji+1 
    429                   pt3d(ji, 1 ,jk) = 0.e0 
     438                  pt3d(ji, 1 ,jk) = zland 
    430439                  pt3d(ji,jpj,jk) = psgn * pt3d(iju,jpj-2,jk) 
    431440               END DO 
    … …  
    437446                  DO ji = 2, jpi 
    438447                     ijt = jpi-ji+2 
    439                      pt3d(ji,  1  ,jk) = 0.e0 
     448                     pt3d(ji,  1  ,jk) = zland 
    440449                     pt3d(ji,jpj-1,jk) = psgn * pt3d(ijt,jpj-2,jk) 
    441450                     pt3d(ji,jpj  ,jk) = psgn * pt3d(ijt,jpj-3,jk) 
    … …  
    451460         CASE ( 5 , 6 )                        ! *  North fold  F-point pivot 
    452461 
    453             pt3d( 1 ,jpj,jk) = 0.e0 
    454             pt3d(jpi,jpj,jk) = 0.e0 
     462            pt3d( 1 ,jpj,jk) = zland 
     463            pt3d(jpi,jpj,jk) = zland 
    455464 
    456465            SELECT CASE ( cd_type ) 
    … …  
    458467               DO ji = 1, jpi 
    459468                  ijt = jpi-ji+1 
    460                   pt3d(ji, 1 ,jk) = 0.e0 
     469                  pt3d(ji, 1 ,jk) = zland 
    461470                  pt3d(ji,jpj,jk) = psgn * pt3d(ijt,jpj-1,jk) 
    462471               END DO 
    … …  
    464473                  DO ji = 1, jpi-1 
    465474                     iju = jpi-ji 
    466                      pt3d(ji, 1 ,jk) = 0.e0 
     475                     pt3d(ji, 1 ,jk) = zland 
    467476                     pt3d(ji,jpj,jk) = psgn * pt3d(iju,jpj-1,jk) 
    468477                  END DO 
    … …  
    470479                  DO ji = 1, jpi 
    471480                     ijt = jpi-ji+1 
    472                      pt3d(ji, 1 ,jk) = 0.e0 
     481                     pt3d(ji, 1 ,jk) = zland 
    473482                     pt3d(ji,jpj,jk) = psgn * pt3d(ijt,jpj-2,jk) 
    474483                  END DO 
    … …  
    492501            SELECT CASE ( cd_type ) 
    493502            CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points 
    494                pt3d(:, 1 ,jk) = 0.e0 
    495                pt3d(:,jpj,jk) = 0.e0 
    496             CASE ( 'F' )                               ! F-point 
    497                pt3d(:,jpj,jk) = 0.e0 
     503               pt3d(:, 1 ,jk) = zland 
     504               pt3d(:,jpj,jk) = zland 
     505            CASE ( 'F' )                               ! F-point 
     506               pt3d(:,jpj,jk) = zland 
    498507            END SELECT 
    499508 
    … …  
    506515 
    507516 
    508    SUBROUTINE lbc_lnk_2d( pt2d, cd_type, psgn, cd_mpp ) 
     517   SUBROUTINE lbc_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval ) 
    509518      !!--------------------------------------------------------------------- 
    510519      !!                 ***  ROUTINE lbc_lnk_2d  *** 
    … …  
    532541      CHARACTER(len=3), INTENT( in ), OPTIONAL ::    & 
    533542         cd_mpp        ! fill the overlap area only (here do nothing) 
     543      REAL(wp)        , INTENT(in   ), OPTIONAL           ::   pval      ! background value (used at closed boundaries) 
    534544 
    535545      !! * Local declarations 
    536546      INTEGER  ::   ji 
    537547      INTEGER  ::   ijt, iju 
     548      REAL(wp) ::   zland 
    538549      !!---------------------------------------------------------------------- 
    539       !!  OPA 8.5, LODYC-IPSL (2002) 
    540       !!---------------------------------------------------------------------- 
     550 
     551      IF( PRESENT( pval ) ) THEN      ! set land value (zero by default) 
     552         zland = pval 
     553      ELSE 
     554         zland = 0.e0 
     555      ENDIF 
    541556 
    542557      IF (PRESENT(cd_mpp)) THEN 
    … …  
    556571         SELECT CASE ( cd_type ) 
    557572         CASE ( 'T' , 'U' , 'V' , 'W' )                ! T-, U-, V-, W-points 
    558             pt2d( 1 ,:) = 0.e0 
    559             pt2d(jpi,:) = 0.e0 
     573            pt2d( 1 ,:) = zland 
     574            pt2d(jpi,:) = zland 
    560575         CASE ( 'F' )                                  ! F-point, ice U-V point 
    561             pt2d(jpi,:) = 0.e0  
     576            pt2d(jpi,:) = zland 
    562577         CASE ( 'I' )                                  ! F-point, ice U-V point 
    563             pt2d( 1 ,:) = 0.e0  
    564             pt2d(jpi,:) = 0.e0  
     578            pt2d( 1 ,:) = zland 
     579            pt2d(jpi,:) = zland 
    565580         END SELECT 
    566581 
    … …  
    576591         CASE ( 'T' , 'U' , 'W' )                      ! T-, U-, W-points 
    577592            pt2d(:, 1 ) = pt2d(:,3) 
    578             pt2d(:,jpj) = 0.e0 
     593            pt2d(:,jpj) = zland 
    579594         CASE ( 'V' , 'F' , 'I' )                      ! V-, F-points, ice U-V point 
    580595            pt2d(:, 1 ) = psgn * pt2d(:,2) 
    581             pt2d(:,jpj) = 0.e0 
     596            pt2d(:,jpj) = zland 
    582597         END SELECT 
    583598 
    584599      CASE ( 3 , 4 )                           ! * North fold  T-point pivot 
    585600 
    586          pt2d( 1 , 1 ) = 0.e0        !!!!!  bug gm ??? !Edmee 
    587          pt2d( 1 ,jpj) = 0.e0 
    588          pt2d(jpi,jpj) = 0.e0 
     601         pt2d( 1 , 1 ) = zland       !!!!!  bug gm ??? !Edmee 
     602         pt2d( 1 ,jpj) = zland 
     603         pt2d(jpi,jpj) = zland 
    589604 
    590605         SELECT CASE ( cd_type ) 
    … …  
    593608            DO ji = 2, jpi 
    594609               ijt = jpi-ji+2 
    595                pt2d(ji, 1 ) = 0.e0 
     610               pt2d(ji, 1 ) = zland 
    596611               pt2d(ji,jpj) = psgn * pt2d(ijt,jpj-2) 
    597612            END DO 
    … …  
    604619            DO ji = 1, jpi-1 
    605620               iju = jpi-ji+1 
    606                pt2d(ji, 1 ) = 0.e0 
     621               pt2d(ji, 1 ) = zland 
    607622               pt2d(ji,jpj) = psgn * pt2d(iju,jpj-2) 
    608623            END DO 
    … …  
    615630            DO ji = 2, jpi 
    616631               ijt = jpi-ji+2 
    617                pt2d(ji, 1   ) = 0.e0 
     632               pt2d(ji, 1   ) = zland 
    618633               pt2d(ji,jpj-1) = psgn * pt2d(ijt,jpj-2) 
    619634               pt2d(ji,jpj  ) = psgn * pt2d(ijt,jpj-3) 
    … …  
    628643 
    629644         CASE ( 'I' )                                  ! ice U-V point 
    630             pt2d(:, 1 ) = 0.e0 
     645            pt2d(:, 1 ) = zland 
    631646            pt2d(2,jpj) = psgn * pt2d(3,jpj-1) 
    632647            DO ji = 3, jpi 
    … …  
    639654      CASE ( 5 , 6 )                           ! * North fold  F-point pivot 
    640655 
    641          pt2d( 1 , 1 ) = 0.e0           !!bug  ??? 
    642          pt2d( 1 ,jpj) = 0.e0 
    643          pt2d(jpi,jpj) = 0.e0 
     656         pt2d( 1 , 1 ) = zland          !!bug  ??? 
     657         pt2d( 1 ,jpj) = zland 
     658         pt2d(jpi,jpj) = zland 
    644659 
    645660         SELECT CASE ( cd_type ) 
    … …  
    648663            DO ji = 1, jpi 
    649664               ijt = jpi-ji+1 
    650                pt2d(ji, 1 ) = 0.e0 
     665               pt2d(ji, 1 ) = zland 
    651666               pt2d(ji,jpj) = psgn * pt2d(ijt,jpj-1) 
    652667            END DO 
    … …  
    655670            DO ji = 1, jpi-1 
    656671               iju = jpi-ji 
    657                pt2d(ji, 1 ) = 0.e0 
     672               pt2d(ji, 1 ) = zland 
    658673               pt2d(ji,jpj) = psgn * pt2d(iju,jpj-1) 
    659674            END DO 
    … …  
    662677            DO ji = 1, jpi 
    663678               ijt = jpi-ji+1 
    664                pt2d(ji, 1 ) = 0.e0 
     679               pt2d(ji, 1 ) = zland 
    665680               pt2d(ji,jpj) = psgn * pt2d(ijt,jpj-2) 
    666681            END DO 
    … …  
    681696 
    682697         CASE ( 'I' )                                  ! ice U-V point 
    683             pt2d( : , 1 ) = 0.e0 
    684             pt2d( 2 ,jpj) = 0.e0 
     698            pt2d( : , 1 ) = zland 
     699            pt2d( 2 ,jpj) = zland 
    685700            DO ji = 2 , jpim1 
    686701               ijt = jpi - ji + 2 
    … …  
    694709         SELECT CASE ( cd_type ) 
    695710         CASE ( 'T' , 'U' , 'V' , 'W' )                ! T-, U-, V-, W-points 
    696             pt2d(:, 1 ) = 0.e0 
    697             pt2d(:,jpj) = 0.e0 
     711            pt2d(:, 1 ) = zland 
     712            pt2d(:,jpj) = zland 
    698713         CASE ( 'F' )                                  ! F-point 
    699             pt2d(:,jpj) = 0.e0 
     714            pt2d(:,jpj) = zland 
    700715         CASE ( 'I' )                                  ! ice U-V point 
    701             pt2d(:, 1 ) = 0.e0 
    702             pt2d(:,jpj) = 0.e0 
     716            pt2d(:, 1 ) = zland 
     717            pt2d(:,jpj) = zland 
    703718         END SELECT 
    704719 

  • trunk/NEMO/OPA_SRC/lib_mpp.F90

    r699 r717  
    597597#endif 
    598598 
    599    SUBROUTINE mpp_lnk_3d( ptab, cd_type, psgn, cd_mpp ) 
     599   SUBROUTINE mpp_lnk_3d( ptab, cd_type, psgn, cd_mpp, pval ) 
    600600      !!---------------------------------------------------------------------- 
    601601      !!                  ***  routine mpp_lnk_3d  *** 
    … …  
    632632      CHARACTER(len=3), INTENT( in ), OPTIONAL ::    & 
    633633         cd_mpp        ! fill the overlap area only  
     634      REAL(wp)        , INTENT(in   ), OPTIONAL           ::   pval      ! background value (used at closed boundaries) 
    634635 
    635636      !! * Local variables 
    … …  
    638639      INTEGER ::   ml_req1, ml_req2, ml_err     ! for key_mpi_isend 
    639640      INTEGER ::   ml_stat(MPI_STATUS_SIZE)     ! for key_mpi_isend 
     641      REAL(wp) ::   zland 
    640642      !!---------------------------------------------------------------------- 
    641643 
    642644      ! 1. standard boundary treatment 
    643645      ! ------------------------------ 
     646 
     647      IF( PRESENT( pval ) ) THEN      ! set land value (zero by default) 
     648         zland = pval 
     649      ELSE 
     650         zland = 0.e0 
     651      ENDIF 
    644652 
    645653      IF( PRESENT( cd_mpp ) ) THEN 
    … …  
    662670            SELECT CASE ( cd_type ) 
    663671            CASE ( 'T', 'U', 'V', 'W' ) 
    664                ptab(     1       :jpreci,:,:) = 0.e0 
    665                ptab(nlci-jpreci+1:jpi   ,:,:) = 0.e0 
     672               ptab(     1       :jpreci,:,:) = zland 
     673               ptab(nlci-jpreci+1:jpi   ,:,:) = zland 
    666674            CASE ( 'F' ) 
    667                ptab(nlci-jpreci+1:jpi   ,:,:) = 0.e0 
     675               ptab(nlci-jpreci+1:jpi   ,:,:) = zland 
    668676            END SELECT  
    669677         ENDIF 
    … …  
    673681         SELECT CASE ( cd_type ) 
    674682         CASE ( 'T', 'U', 'V', 'W' ) 
    675             ptab(:,     1       :jprecj,:) = 0.e0 
    676             ptab(:,nlcj-jprecj+1:jpj   ,:) = 0.e0 
     683            ptab(:,     1       :jprecj,:) = zland 
     684            ptab(:,nlcj-jprecj+1:jpj   ,:) = zland 
    677685         CASE ( 'F' ) 
    678             ptab(:,nlcj-jprecj+1:jpj   ,:) = 0.e0 
     686            ptab(:,nlcj-jprecj+1:jpj   ,:) = zland 
    679687         END SELECT 
    680688      
    … …  
    10501058 
    10511059 
    1052    SUBROUTINE mpp_lnk_2d( pt2d, cd_type, psgn, cd_mpp ) 
     1060   SUBROUTINE mpp_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval ) 
    10531061      !!---------------------------------------------------------------------- 
    10541062      !!                  ***  routine mpp_lnk_2d  *** 
    … …  
    10841092      CHARACTER(len=3), INTENT( in ), OPTIONAL ::    & 
    10851093         cd_mpp        ! fill the overlap area only  
     1094      REAL(wp)        , INTENT(in   ), OPTIONAL           ::   pval      ! background value (used at closed boundaries) 
    10861095 
    10871096      !! * Local variables 
    … …  
    10921101      INTEGER  ::   ml_req1, ml_req2, ml_err     ! for key_mpi_isend 
    10931102      INTEGER  ::   ml_stat(MPI_STATUS_SIZE)     ! for key_mpi_isend 
     1103      REAL(wp) ::   zland 
    10941104      !!---------------------------------------------------------------------- 
     1105 
     1106      IF( PRESENT( pval ) ) THEN      ! set land value (zero by default) 
     1107         zland = pval 
     1108      ELSE 
     1109         zland = 0.e0 
     1110      ENDIF 
    10951111 
    10961112      ! 1. standard boundary treatment 
    … …  
    11151131            SELECT CASE ( cd_type ) 
    11161132            CASE ( 'T', 'U', 'V', 'W' , 'I' ) 
    1117                pt2d(     1       :jpreci,:) = 0.e0 
    1118                pt2d(nlci-jpreci+1:jpi   ,:) = 0.e0 
     1133               pt2d(     1       :jpreci,:) = zland 
     1134               pt2d(nlci-jpreci+1:jpi   ,:) = zland 
    11191135            CASE ( 'F' ) 
    1120                pt2d(nlci-jpreci+1:jpi   ,:) = 0.e0 
     1136               pt2d(nlci-jpreci+1:jpi   ,:) = zland 
    11211137            END SELECT 
    11221138         ENDIF 
    … …  
    11261142         SELECT CASE ( cd_type ) 
    11271143         CASE ( 'T', 'U', 'V', 'W' , 'I' ) 
    1128             pt2d(:,     1       :jprecj) = 0.e0 
    1129             pt2d(:,nlcj-jprecj+1:jpj   ) = 0.e0 
     1144            pt2d(:,     1       :jprecj) = zland 
     1145            pt2d(:,nlcj-jprecj+1:jpj   ) = zland 
    11301146         CASE ( 'F' ) 
    1131             pt2d(:,nlcj-jprecj+1:jpj   ) = 0.e0 
     1147            pt2d(:,nlcj-jprecj+1:jpj   ) = zland 
    11321148         END SELECT 
    11331149 
    … …  
    13941410   
    13951411            CASE ( 'I' )                                  ! ice U-V point 
    1396                pt2d( 2 ,nlcj) = 0.e0 
     1412               pt2d( 2 ,nlcj) = zland 
    13971413               DO ji = 2 , nlci-1 
    13981414                  ijt = iloc - ji + 2 

  • trunk/NEMO/OPA_SRC/opa.F90

    r709 r717  
    309309#endif 
    310310 
    311       CALL flx_init                         ! Thermohaline forcing initialization 
    312  
    313       CALL flx_fwb_init                     ! FreshWater Budget correction 
    314  
    315311      CALL dia_ptr_init                     ! Poleward TRansports initialization 
    316312 

  • trunk/NEMO/OPA_SRC/restart.F90

    r709 r717  
    1919   USE phycst          ! physical constants 
    2020   USE daymod          ! calendar 
    21    USE ice_oce         ! ice variables 
    2221   USE cpl_oce, ONLY : lk_cpl              ! 
    2322   USE in_out_manager  ! I/O manager 
    … …  
    138137      CALL iom_rstput( kt, nitrst, numrow, 'rotn'   , rotn    ) 
    139138      CALL iom_rstput( kt, nitrst, numrow, 'hdivn'  , hdivn   ) 
    140  
    141 #if defined key_ice_lim         
    142       CALL iom_rstput( kt, nitrst, numrow, 'nfice'  , REAL( nfice, wp) )   !  ice computation frequency 
    143       CALL iom_rstput( kt, nitrst, numrow, 'sst_io' , sst_io  ) 
    144       CALL iom_rstput( kt, nitrst, numrow, 'sss_io' , sss_io  ) 
    145       CALL iom_rstput( kt, nitrst, numrow, 'u_io'   , u_io    ) 
    146       CALL iom_rstput( kt, nitrst, numrow, 'v_io'   , v_io    ) 
    147 # if defined key_coupled 
    148       CALL iom_rstput( kt, nitrst, numrow, 'alb_ice', alb_ice ) 
    149 # endif 
    150 #endif 
    151 #if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core  
    152       CALL iom_rstput( kt, nitrst, numrow, 'nfbulk' , REAL( nfbulk, wp) )   !  bulk computation frequency 
    153       CALL iom_rstput( kt, nitrst, numrow, 'gsst'   , gsst    ) 
    154 #endif 
    155139 
    156140      IF( nn_dynhpg_rst == 1 .OR. lk_vvl ) THEN 
    … …  
    204188      !!                    has been stored in the restart file. 
    205189      !!---------------------------------------------------------------------- 
    206       REAL(wp) ::   zcoef, zkt, zrdt, zrdttra1, zndastp, znfice, znfbulk 
     190      REAL(wp) ::   zcoef, zkt, zrdt, zrdttra1, zndastp 
    207191#if defined key_ice_lim 
    208192      INTEGER  ::   ji, jj 
    … …  
    299283      ENDIF 
    300284 
    301       !!sm: TO BE MOVED IN NEW SURFACE MODULE... 
    302  
    303 #if defined key_ice_lim 
    304       ! Louvain La Neuve Sea Ice Model 
    305       IF( iom_varid( numror, 'nfice' ) > 0 ) then  
    306          CALL iom_get( numror             , 'nfice'  , znfice  )   ! ice computation frequency 
    307          CALL iom_get( numror, jpdom_autoglo, 'sst_io' , sst_io  ) 
    308          CALL iom_get( numror, jpdom_autoglo, 'sss_io' , sss_io  ) 
    309          CALL iom_get( numror, jpdom_autoglo, 'u_io'   , u_io    ) 
    310          CALL iom_get( numror, jpdom_autoglo, 'v_io'   , v_io    ) 
    311 # if defined key_coupled 
    312          CALL iom_get( numror, jpdom_autoglo, 'alb_ice', alb_ice ) 
    313 # endif 
    314          IF( znfice /= REAL( nfice, wp ) ) THEN      ! if nfice changed between 2 runs 
    315             zcoef = REAL( nfice-1, wp ) / znfice 
    316             sst_io(:,:) = zcoef * sst_io(:,:) 
    317             sss_io(:,:) = zcoef * sss_io(:,:) 
    318             u_io  (:,:) = zcoef * u_io  (:,:) 
    319             v_io  (:,:) = zcoef * v_io  (:,:) 
    320          ENDIF 
    321       ELSE 
    322          IF(lwp) WRITE(numout,*) 
    323          IF(lwp) WRITE(numout,*) 'rst_read :  LLN sea Ice Model => Ice initialization' 
    324          IF(lwp) WRITE(numout,*) 
    325          zcoef = REAL( nfice-1, wp ) 
    326          sst_io(:,:) = zcoef *( tn(:,:,1) + rt0 )          !!bug a explanation is needed here! 
    327          sss_io(:,:) = zcoef *  sn(:,:,1) 
    328          zcoef = 0.5 * REAL( nfice-1, wp ) 
    329          DO jj = 2, jpj 
    330             DO ji = fs_2, jpi   ! vector opt. 
    331                u_io(ji,jj) = zcoef * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) ) 
    332                v_io(ji,jj) = zcoef * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) ) 
    333             END DO 
    334          END DO 
    335 # if defined key_coupled 
    336          alb_ice(:,:) = 0.8 * tmask(:,:,1) 
    337 # endif 
    338       ENDIF 
    339 #endif 
    340 #if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core  
    341       ! Louvain La Neuve Sea Ice Model 
    342       IF( iom_varid( numror, 'nfbulk' ) > 0 ) THEN  
    343          CALL iom_get( numror             , 'nfbulk', znfbulk )   ! bulk computation frequency 
    344          CALL iom_get( numror, jpdom_autoglo, 'gsst'  , gsst    ) 
    345          IF( znfbulk /= REAL(nfbulk, wp) ) THEN      ! if you change nfbulk between 2 runs 
    346             zcoef = REAL( nfbulk-1, wp ) / znfbulk 
    347             gsst(:,:) = zcoef * gsst(:,:) 
    348          ENDIF 
    349       ELSE 
    350          IF(lwp) WRITE(numout,*) 
    351          IF(lwp) WRITE(numout,*) 'rst_read :  LLN sea Ice Model => Ice initialization' 
    352          IF(lwp) WRITE(numout,*) 
    353          gsst(:,:) = REAL( nfbulk - 1, wp )*( tn(:,:,1) + rt0 ) 
    354       ENDIF 
    355 #endif 
    356  
    357       !!sm: end of TO BE MOVED IN NEW SURFACE MODULE... 
    358  
    359285      IF( iom_varid( numror, 'rhd' ) > 0 ) THEN 
    360286         CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd  ) 
Note: See TracChangeset for help on using the changeset viewer.