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

Changeset 1857


Ignore:
Timestamp:
2010-05-03T13:59:46+02:00 (14 years ago)
Author:
gm
Message:

ticket:#665 Reverting previous commit and going back to revision 1850

Location:
branches/DEV_r1837_mass_heat_salt_fluxes/NEMO
Files:
21 edited
1 copied

Legend:

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

    r1855 r1857  
    22   !!====================================================================== 
    33   !!                   ***  MODULE  dom_ice  *** 
    4    !! LIM-2 Sea Ice :   Domain  variables 
     4   !! LIM 2.0 Sea Ice :   Domain  variables 
    55   !!====================================================================== 
    6    !! History :   2.0  !  2003-08  (C. Ethe)  Free form and module 
     6   !! History :   2.0  !  03-08  (C. Ethe)  Free form and module 
    77   !!---------------------------------------------------------------------- 
    88#if defined key_lim2 
    99   !!---------------------------------------------------------------------- 
    10    !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
     10   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005) 
     11   !! $Id$ 
     12   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    1113   !!---------------------------------------------------------------------- 
    1214   USE par_ice_2 
     
    1618 
    1719   LOGICAL, PUBLIC ::   l_jeq     = .TRUE.     !: Equator inside the domain flag 
     20 
    1821   INTEGER, PUBLIC ::   njeq , njeqm1          !: j-index of the equator if it is inside the domain 
     22      !                                        !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
    1923 
    20    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor, fcor     !: coriolis factor and coeficient 
    21    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   covrai           !: sine of geographic latitude 
    22    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   area             !: surface of grid cell  
    23    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tms   , tmu      !: temperature and velocity points masks 
    24    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght             !: weight of the 4 neighbours to compute averages 
    25    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   akappa, bkappa   !: first and third group of metric coefficients 
    26    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd           !: second group of metric coefficients 
     24   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor,   &  !: coriolis factor and coeficient 
     25      &                                              covrai ,         &  !: sine of geographic latitude 
     26      &                                              area   ,         &  !: surface of grid cell  
     27      &                                              tms    , tmu        !: temperature and velocity points masks 
     28   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght   ,         &  !: weight of the 4 neighbours to compute averages 
     29      &                                              akappa , bkappa     !: first and third group of metric coefficients 
     30   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd   !: second group of metric coefficients 
    2731 
    28 #else 
    29    !!---------------------------------------------------------------------- 
    30    !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    31    !!---------------------------------------------------------------------- 
     32   !!====================================================================== 
    3233#endif 
    33  
    34    !!-------- ------------------------------------------------------------- 
    35    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
    36    !! $Id$ 
    37    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    38    !!====================================================================== 
    3934END MODULE dom_ice_2 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90

    r1855 r1857  
    44   !! Sea Ice physics:  diagnostics variables of ice defined in memory 
    55   !!===================================================================== 
    6    !! History :  2.0  !  2003-08 (C. Ethe)  F90: Free form and module 
    7    !!             -   !  2010-05 (G. Madec) suppression of thd_ice_2 
     6   !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module 
    87   !!---------------------------------------------------------------------- 
    98#if defined key_lim2 
     
    1110   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1211   !!---------------------------------------------------------------------- 
     12   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006) 
     13   !! $Id$ 
     14   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     15   !!---------------------------------------------------------------------- 
     16   !! * Modules used 
    1317   USE par_ice_2          ! LIM sea-ice parameters 
    1418 
     
    1620   PRIVATE 
    1721 
    18    !---------------------------------------------------! 
    19    !                 Sea-Ice namelist                  ! 
    20    !---------------------------------------------------! 
    21    !                                                                !!! **  Namelist namicerun  ** 
    22    CHARACTER(len=32), PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input) 
    23    CHARACTER(len=32), PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output) 
    24    LOGICAL          , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F) 
    25    LOGICAL          , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping 
    26    REAL(wp)         , PUBLIC ::   hsndif        = 0.e0               !: computation of temp. in snow (0) or not (9999) 
    27    REAL(wp)         , PUBLIC ::   hicdif        = 0.e0               !: computation of temp. in ice (0) or not (9999) 
    28    REAL(wp)         , PUBLIC ::   acrit(2) = (/ 1.e-06 , 1.e-06 /)   !: min lead fraction in north and south hemispheres 
    29     
    30    !                                         !!! **  Namelist  namicedyn  ** 
    31    INTEGER , PUBLIC ::   nbiter =   1         !: number of sub-time steps for relaxation 
    32    INTEGER , PUBLIC ::   nbitdr = 250         !: maximum number of iterations for relaxation 
    33    REAL(wp), PUBLIC ::   rdt_ice              !: ice time step 
    34    REAL(wp), PUBLIC ::   epsd   =   1.0e-20   !: tolerance parameter for dynamic 
    35    REAL(wp), PUBLIC ::   alpha  =   0.5       !: coefficient for semi-implicit coriolis 
    36    REAL(wp), PUBLIC ::   dm     =   0.6e+03   !: diffusion constant for dynamics 
    37    REAL(wp), PUBLIC ::   om     =   0.5       !: relaxation constant 
    38    REAL(wp), PUBLIC ::   resl   =   5.0e-05   !: maximum value for the residual of relaxation 
    39    REAL(wp), PUBLIC ::   cw     =   5.0e-03   !: drag coefficient for oceanic stress 
    40    REAL(wp), PUBLIC ::   angvg  =   0.e0      !: turning angle for oceanic stress 
    41    REAL(wp), PUBLIC ::   pstar  =   1.0e+04   !: first bulk-rheology parameter 
    42    REAL(wp), PUBLIC ::   c_rhg  =  20.e0      !: second bulk-rhelogy parameter 
    43    REAL(wp), PUBLIC ::   etamn  =   0.e+07    !: minimun value for viscosity 
    44    REAL(wp), PUBLIC ::   creepl =   2.e-08    !: creep limit 
    45    REAL(wp), PUBLIC ::   ecc    =   2.e0      !: eccentricity of the elliptical yield curve 
    46    REAL(wp), PUBLIC ::   ahi0   = 350.e0      !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
     22   !!* Share parameters namelist (namicerun read in iceini) * 
     23   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input) 
     24   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output) 
     25   LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F) 
     26   LOGICAL               , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping 
     27   REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of temp. in snow (0) or not (9999) 
     28   REAL(wp)              , PUBLIC ::   hicdif        = 0.e0               !: computation of temp. in ice (0) or not (9999) 
     29   REAL(wp), DIMENSION(2), PUBLIC ::   acrit = (/ 1.e-06 , 1.e-06 /)      !: minimum fraction for leads in  
     30   !                                                                      !: north and south hemisphere 
     31   !!* ice-dynamic namelist (namicedyn) * 
     32   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation 
     33   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation 
     34   REAL(wp), PUBLIC ::   rdt_ice            !: ice time step 
     35   REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic 
     36   REAL(wp), PUBLIC ::   alpha  = 0.5       !: coefficient for semi-implicit coriolis 
     37   REAL(wp), PUBLIC ::   dm     = 0.6e+03   !: diffusion constant for dynamics 
     38   REAL(wp), PUBLIC ::   om     = 0.5       !: relaxation constant 
     39   REAL(wp), PUBLIC ::   resl   = 5.0e-05   !: maximum value for the residual of relaxation 
     40   REAL(wp), PUBLIC ::   cw     = 5.0e-03   !: drag coefficient for oceanic stress 
     41   REAL(wp), PUBLIC ::   angvg  = 0.e0      !: turning angle for oceanic stress 
     42   REAL(wp), PUBLIC ::   pstar  = 1.0e+04   !: first bulk-rheology parameter 
     43   REAL(wp), PUBLIC ::   c_rhg  = 20.e0     !: second bulk-rhelogy parameter 
     44   REAL(wp), PUBLIC ::   etamn  = 0.e+07    !: minimun value for viscosity 
     45   REAL(wp), PUBLIC ::   creepl = 2.e-08    !: creep limit 
     46   REAL(wp), PUBLIC ::   ecc    = 2.e0      !: eccentricity of the elliptical yield curve 
     47   REAL(wp), PUBLIC ::   ahi0   = 350.e0    !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    4748 
    48    !                                              !!! ** Namelist  namicethd  ** 
    49    REAL(wp), PUBLIC ::   hmelt      = -0.15        !: maximum melting at the bottom 
    50    REAL(wp), PUBLIC ::   hicmin     =  0.2         !: ice th. corr. to max. ener. in brine pocket 
    51    REAL(wp), PUBLIC ::   hiclim     =  0.05        !: minimum ice thickness 
    52    REAL(wp), PUBLIC ::   amax       =  0.999       !: maximum lead fraction 
    53    REAL(wp), PUBLIC ::   swiqst     =  1.0         !: energy stored in brine pocket (1) or not (0) 
    54    REAL(wp), PUBLIC ::   sbeta      =  1.0         !: numerical scheme for diffusion in ice  
    55    REAL(wp), PUBLIC ::   parlat     =  0.0         !: percent. of energy used for lateral ablation 
    56    REAL(wp), PUBLIC ::   hakspl     =  0.5         !: slope of distr. for Hakkinen-Mellro's lat. melt 
    57    REAL(wp), PUBLIC ::   hibspl     =  0.5         !: slope of distribution for Hibler's lat. melt 
    58    REAL(wp), PUBLIC ::   exld       =  2.0         !: exponent for leads-closure rate 
    59    REAL(wp), PUBLIC ::   hakdif     =  1.0         !: coefficient for diffusions of ice and snow 
    60    REAL(wp), PUBLIC ::   thth       =  0.2         !: thick. for comp. of eq. thermal conduct 
    61    REAL(wp), PUBLIC ::   hnzst      =  0.1         !: thick. of the surf. layer in temp. comp. 
    62    REAL(wp), PUBLIC ::   parsub     =  1.0         !: switch for snow sublimation or not 
    63    REAL(wp), PUBLIC ::   alphs      =  1.0         !: coef. for snow density when snow-ice formation 
    64    REAL(wp), PUBLIC ::   hiccrit(2) = (/0.3,0.3/)  !: ice th. for lateral accretion in the NH (SH) (m) 
    65  
    66    !---------------------------------------------------! 
    67    !                Sea-Ice variables                  ! 
    68    !---------------------------------------------------! 
    69  
    70    REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( ecc * ecc ) 
    71    REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cw 
    72    REAL(wp), PUBLIC ::   sangvg, cangvg   !: sin and cos of the turning angle for ocean stress 
    73    REAL(wp), PUBLIC ::   pstarh           !: pstar / 2.0 
    74    REAL(wp), PUBLIC ::   uscomi           !: inverse of minimum lead fraction      !!gm to be suppress 
    75    REAL(wp), PUBLIC ::   cnscg            !: ratio  rcpsn/rcpic                    !!gm to be suppress 
     49   REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc ) 
     50   REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw 
     51   REAL(wp), PUBLIC ::   sangvg, cangvg     !: sin and cos of the turning angle for ocean stress 
     52   REAL(wp), PUBLIC ::   pstarh             !: pstar / 2.0 
    7653 
    7754   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
    7855   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
    7956   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses 
    80    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s         !: friction velocity 
     57   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity 
    8158 
    82    INTEGER , PUBLIC, DIMENSION(jpij) ::   npb         !: number of points where computations has to be done 
    83    INTEGER , PUBLIC, DIMENSION(jpij) ::   npac        !: correspondance between the points 
     59   !!* diagnostic quantities 
     60   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin) 
     61   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS 
     62   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness 
     63   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnif         !: Snow thickness 
     64   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicifp        !: Ice production/melting 
     65   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   frld          !: Leads fraction = 1-a/totalarea 
     66   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   phicif        !: ice thickness  at previous time  
     67   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pfrld         !: Leads fraction at previous time   
     68   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qstoif        !: Energy stored in the brine pockets 
     69   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fbif          !: Heat flux at the ice base 
     70   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmsnif       !: Variation of snow mass 
     71   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass 
     72   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean) 
     73   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing  
     74   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean 
     75   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2) 
     76   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget 
     77   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice 
     78   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max heat contained in brine pockets (?) 
     79   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?) 
     80   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?) 
     81   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy in case of toral lateral ablation (?) 
     82   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice 
    8483 
    85    !                                                         !!* sea-ice variables * 
    86    REAL(wp), PUBLIC ::   sist   (jpi,jpj), sist_1d   (jpij)   !: Sea-Ice Surface Temperature (Kelvin) 
    87    REAL(wp), PUBLIC ::   tfu    (jpi,jpj), tfu_1d    (jpij)   !: Freezing/Melting point temperature of sea water at SSS 
    88    REAL(wp), PUBLIC ::   hicif  (jpi,jpj), h_ice_1d  (jpij)   !: ice thickness 
    89    REAL(wp), PUBLIC ::   phicif (jpi,jpj)                     !: ice thickness  at previous time  
    90    REAL(wp), PUBLIC ::   hsnif  (jpi,jpj), h_snow_1d (jpij)   !: Snow thickness 
    91    REAL(wp), PUBLIC ::   hicifp (jpi,jpj)                     !: Ice production/melting 
    92    REAL(wp), PUBLIC ::   frld   (jpi,jpj), frld_1d   (jpij)   !: Leads fraction = 1-a/totalarea 
    93    REAL(wp), PUBLIC ::   pfrld  (jpi,jpj)                     !: Leads fraction at previous time   
    94    REAL(wp), PUBLIC ::   qstoif (jpi,jpj), qstbif_1d (jpij)   !: Energy stored in the brine pockets         !!gm  
    95    REAL(wp), PUBLIC ::   fbif   (jpi,jpj), fbif_1d   (jpij)   !: Heat flux at the ice base 
    96    REAL(wp), PUBLIC ::   rdmsnif(jpi,jpj), rdmsnif_1d(jpij)   !: Variation of snow mass 
    97    REAL(wp), PUBLIC ::   rdmicif(jpi,jpj), rdmicif_1d(jpij)   !: Variation of ice mass 
    98    REAL(wp), PUBLIC ::   qldif  (jpi,jpj), qldif_1d  (jpij)   !: heat balance of the lead (or of the open ocean) 
    99    REAL(wp), PUBLIC ::   qcmif  (jpi,jpj), qcmif_1d  (jpij)   !: Energy needed to freeze the ocean surface layer  
    100    REAL(wp), PUBLIC ::   fdtcn  (jpi,jpj)                     !: net downward heat flux from the ice to the ocean 
    101    REAL(wp), PUBLIC ::   qdtcn  (jpi,jpj)                     !: energy from the ice to the ocean point (at a factor 2) 
    102    REAL(wp), PUBLIC ::   thcm   (jpi,jpj), thcm_1d   (jpij)   !: part of the solar energy used in the lead heat budget 
    103    REAL(wp), PUBLIC ::   fstric (jpi,jpj), fstric_1d (jpij)   !: Solar flux transmitted trough the ice 
    104    REAL(wp), PUBLIC ::   ffltbif(jpi,jpj), fltbif_1d (jpij)   !: Array linked with the max heat contained in brine pockets (?) 
    105    REAL(wp), PUBLIC ::   fscmbq (jpi,jpj), fscbq_1d  (jpij)   !: Linked with the solar flux below the ice (?) 
    106    REAL(wp), PUBLIC ::   fsbbq  (jpi,jpj)                     !: Also linked with the solar flux below the ice (?) 
    107    REAL(wp), PUBLIC ::   qfvbq  (jpi,jpj), qfvbq_1d  (jpij)   !: Array used to store energy in case of total lateral ablation (?) 
    108    REAL(wp), PUBLIC ::   dmgwi(jpi,jpj), dmgwi_1d(jpij)       !: Variation of the mass of snow ice 
     84   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice   !: two components of the ice   velocity at I-point (m/s) 
     85   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce   !: two components of the ocean velocity at I-point (m/s) 
    10986 
    110    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice     !: two components of the ice   velocity at I-point (m/s) 
    111    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce     !: two components of the ocean velocity at I-point (m/s) 
     87   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer 
    11288 
    113    REAL(wp), PUBLIC ::   tbif(jpi,jpj,jplayersp1), tbif_1d(jpij,jplayersp1)  !: Temperature inside the ice+snow layer 
    114  
    115    ! Surface forcing transford into  1D field  (2D field defined in /OPA_SRC/SBC/sbc_ice.F90) 
    116    REAL(wp), PUBLIC, DIMENSION(jpij) ::   qns_ice_1d   !:   qns_ice    
    117    REAL(wp), PUBLIC, DIMENSION(jpij) ::   qsr_ice_1d   !:   qsr_ice 
    118    REAL(wp), PUBLIC, DIMENSION(jpij) ::   qla_ice_1d   !:   qla_ice 
    119    REAL(wp), PUBLIC, DIMENSION(jpij) ::   dqla_ice_1d  !:   dqla_ice 
    120    REAL(wp), PUBLIC, DIMENSION(jpij) ::   dqns_ice_1d  !:   dqns_ice 
    121    REAL(wp), PUBLIC, DIMENSION(jpij) ::   fr1_i0_1d    !:   fr1_i0 
    122    REAL(wp), PUBLIC, DIMENSION(jpij) ::   fr2_i0_1d    !:   fr2_i0 
    123  
    124    ! Surface forcing transformed into 1D field  (2D field defined in /OPA_SRC/SBC/sbc_oce.F90) 
    125    REAL(wp), PUBLIC, DIMENSION(jpij) ::   sprecip_1d   !:   sprecip 
    126     
    127    ! local variable  transformed into 1D field  (2D field defined in /LIM_SRC_2/limthd_2.F90) 
    128    REAL(wp), PUBLIC, DIMENSION(jpij) ::   qlbbq_1d     !:   zqlbsbq 
    129    REAL(wp), PUBLIC, DIMENSION(jpij) ::   dvsbq_1d     !:   zdvosif 
    130    REAL(wp), PUBLIC, DIMENSION(jpij) ::   dvbbq_1d     !:   zdvobif 
    131    REAL(wp), PUBLIC, DIMENSION(jpij) ::   dvlbq_1d     !:   zdvolif 
    132    REAL(wp), PUBLIC, DIMENSION(jpij) ::   dvnbq_1d     !:   zdvonif 
    133     
    134     
    135    REAL(wp), PUBLIC, DIMENSION(jpij) ::   rdvomif_1d    !:   rdvomif       <<<== !gm should be suppressed ! 
    136  
    137    !---------------------------------------------------! 
    138    !   Prather' advection scheme : 1st & 2nd moments   ! 
    139    !---------------------------------------------------! 
     89   !!* moment used in the advection scheme 
    14090   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume 
    14191   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume 
     
    152102#endif 
    153103 
    154    !!-------- ------------------------------------------------------------- 
    155    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
    156    !! $Id$ 
    157    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    158104   !!====================================================================== 
    159105END MODULE ice_2 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/iceini_2.F90

    r1855 r1857  
    44   !!   Sea-ice model : LIM 2.0 Sea ice model Initialization 
    55   !!====================================================================== 
    6    !! History :  1.0  !  2002-08  (G. Madec)  F90: Free form and modules 
    7    !!            2.0  !  2003-08  (C. Ethe)  add ice_run 
     6   !! History :   1.0  !  02-08  (G. Madec)  F90: Free form and modules 
     7   !!             2.0  !  03-08  (C. Ethe)  add ice_run 
    88   !!---------------------------------------------------------------------- 
    99#if defined key_lim2 
     
    1111   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
    13    !!   ice_init_2      : sea-ice model initialization 
    1413   !!---------------------------------------------------------------------- 
    15    USE dom_oce          ! ocean domain 
    16    USE dom_ice_2        ! LIM-2 ice domain 
    17    USE sbc_oce          ! surface boundary condition: ocean 
    18    USE sbc_ice          ! surface boundary condition: sea ice 
    19    USE phycst           ! physical constant 
    20    USE ice_2            ! LIM-2 ice varibles 
    21    USE limmsh_2         ! LIM-2 mesh 
    22    USE limistate_2      ! LIM-2 initial state 
    23    USE limrst_2         ! LIM-2 restart 
    24    USE in_out_manager   ! I/O manager 
     14   !!   ice_init_2       : sea-ice model initialization 
     15   !!   ice_run_2        : Definition some run parameter for ice model 
     16   !!---------------------------------------------------------------------- 
     17   USE dom_oce 
     18   USE dom_ice_2 
     19   USE sbc_oce         ! surface boundary condition: ocean 
     20   USE sbc_ice         ! surface boundary condition: ice 
     21   USE phycst          ! Define parameters for the routines 
     22   USE ice_2 
     23   USE limmsh_2 
     24   USE limistate_2 
     25   USE limrst_2    
     26   USE in_out_manager 
    2527       
    2628   IMPLICIT NONE 
    2729   PRIVATE 
    2830 
    29    PUBLIC   ice_init_2   ! called by sbcice_lim_2.F90 
     31   PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90 
    3032 
    3133   !!---------------------------------------------------------------------- 
    32    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     34   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    3335   !! $Id$  
    34    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     36   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3537   !!---------------------------------------------------------------------- 
    3638 
     
    4345      !! ** purpose :    
    4446      !!---------------------------------------------------------------------- 
     47      ! 
     48      ! Open the namelist file  
     49      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )       
     50      CALL ice_run_2                    !  read in namelist some run parameters 
     51                  
     52      ! Louvain la Neuve Ice model 
     53      rdt_ice = nn_fsbc * rdttra(1) 
     54 
     55      CALL lim_msh_2                  ! ice mesh initialization 
     56      
     57      ! Initial sea-ice state 
     58      IF( .NOT.ln_rstart ) THEN 
     59         CALL lim_istate_2            ! start from rest: sea-ice deduced from sst 
     60      ELSE 
     61         CALL lim_rst_read_2          ! start from a restart file 
     62      ENDIF 
     63       
     64      tn_ice(:,:,1) = sist(:,:)         ! initialisation of ice temperature    
     65      fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction     
     66      ! 
     67   END SUBROUTINE ice_init_2 
     68 
     69 
     70   SUBROUTINE ice_run_2 
     71      !!------------------------------------------------------------------- 
     72      !!                  ***  ROUTINE ice_run_2 *** 
     73      !!                  
     74      !! ** Purpose :   Definition some run parameter for ice model 
     75      !! 
     76      !! ** Method  :   Read the namicerun namelist and check the parameter  
     77      !!       values called at the first timestep (nit000) 
     78      !! 
     79      !! ** input   :   Namelist namicerun 
     80      !!------------------------------------------------------------------- 
    4581      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, ln_limdmp, acrit, hsndif, hicdif 
    4682      !!------------------------------------------------------------------- 
    47       ! 
    48       !                                      ! Open the ice namelist file  
    49       CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )     
    50          
    51       REWIND ( numnam_ice )                  ! Read namicerun namelist 
     83      !                     
     84      REWIND ( numnam_ice )                       ! Read Namelist namicerun  
    5285      READ   ( numnam_ice , namicerun ) 
    53       ! 
    54       IF(lwp) THEN                           ! control print 
     86 
     87      IF(lwp) THEN 
    5588         WRITE(numout,*) 
    5689         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
     
    6295         WRITE(numout,*) '   computation of temp. in ice  (=0) or not (=9999) hicdif = ', hicdif 
    6396      ENDIF 
    64                         
    65       rdt_ice = nn_fsbc * rdttra(1)          ! ice time step  
    66  
    67       CALL lim_msh_2                         ! ice mesh initialization 
    68       
    69       !                                      ! Initial sea-ice state 
    70       IF( .NOT.ln_rstart ) THEN   ;   CALL lim_istate_2        ! start from rest: sea-ice deduced from sst 
    71       ELSE                        ;   CALL lim_rst_read_2      ! restart from a file 
    72       ENDIF 
    73  
    74       !                                      ! ice variables see by the ocean       
    75       tn_ice(:,:,1) = sist(:,:)                 ! surface ice temperature    
    76       fr_i  (:,:)   = 1.0 - frld(:,:)           ! sea-ice fraction     
    7797      ! 
    78    END SUBROUTINE ice_init_2 
     98   END SUBROUTINE ice_run_2 
    7999 
    80100#else 
     
    82102   !!   Default option :        Empty module       NO LIM 2.0 sea-ice model 
    83103   !!---------------------------------------------------------------------- 
     104CONTAINS 
     105   SUBROUTINE ice_init_2      ! Empty routine 
     106   END SUBROUTINE ice_init_2 
    84107#endif 
    85108 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limadv_2.F90

    r1855 r1857  
    22   !!====================================================================== 
    33   !!                       ***  MODULE limadv_2   *** 
    4    !! LIM-2 sea-ice model : sea-ice advection 
     4   !! LIM 2.0 sea-ice model : sea-ice advection 
    55   !!====================================================================== 
    66   !! History :  OPA  ! 2000-01 (LIM)  Original code 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdia_2.F90

    r1855 r1857  
    1212   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
     14   !!---------------------------------------------------------------------- 
    1415   !!   lim_dia_2      : computation of the time evolution of keys var. 
    1516   !!   lim_dia_init_2 : initialization and namelist read 
     
    2728   PRIVATE 
    2829 
    29    PUBLIC   lim_dia_2   ! called by sbc_ice_lim_2 
    30     
    31    INTEGER, PUBLIC ::   ntmoy   = 1   !: instantaneous values of ice evolution or averaging ntmoy 
    32    INTEGER, PUBLIC ::   ninfo   = 1   !: frequency of ouputs on file ice_evolu in case of averaging 
    33  
    34    !                                              !!! Parameters for outputs to files "evolu" 
    35    INTEGER, PARAMETER ::   jpinfmx = 100           ! maximum number of key variables 
    36    INTEGER, PARAMETER ::   jpchinf = 5             ! ??? 
    37    INTEGER, PARAMETER ::   jpchsep = jpchinf + 2   ! ??? 
    38  
    39    INTEGER ::   nfrinf  = 4   ! number of variables written in one line  
    40    INTEGER ::   nferme        ! last time step at which the var. are written on file 
    41    INTEGER ::   nvinfo        ! number of total variables  
    42    INTEGER ::   nbvt          ! number of time variables 
    43    INTEGER ::   naveg         ! number of step for accumulation before averaging 
    44  
    45    CHARACTER(len= 8) ::   fmtinf  = '1PE13.5 '   ! format of the output values   
    46    CHARACTER(len=30) ::   fmtw, fmtr, fmtitr     ! various formats 
    47    CHARACTER(len=jpchsep), DIMENSION(jpinfmx) ::   titvar   ! title of key variables 
     30   PUBLIC               lim_dia_2          ! called by sbc_ice_lim_2 
     31   INTEGER, PUBLIC ::   ntmoy   = 1 ,   &  !: instantaneous values of ice evolution or averaging ntmoy 
     32      &                 ninfo   = 1        !: frequency of ouputs on file ice_evolu in case of averaging 
     33 
     34   INTEGER, PARAMETER ::   &  ! Parameters for outputs to files "evolu" 
     35      jpinfmx = 100         ,    &  ! maximum number of key variables 
     36      jpchinf = 5           ,    &  ! ??? 
     37      jpchsep = jpchinf + 2         ! ??? 
     38 
     39   INTEGER ::   & 
     40      nfrinf  = 4 ,     &  ! number of variables written in one line  
     41      nferme ,          &  ! last time step at which the var. are written on file 
     42      nvinfo ,          &  ! number of total variables  
     43      nbvt   ,          &  ! number of time variables 
     44      naveg                ! number of step for accumulation before averaging 
     45 
     46   CHARACTER(len= 8) ::   fmtinf  = '1PE13.5 ' ! format of the output values   
     47   CHARACTER(len=30) ::   fmtw  ,           &  ! formats 
     48      &                   fmtr  ,           &  ! ??? 
     49      &                   fmtitr               ! ??? 
     50   CHARACTER(len=jpchsep), DIMENSION(jpinfmx) ::   titvar               ! title of key variables 
    4851  
    4952   REAL(wp)                     ::   epsi06 = 1.e-06      ! ??? 
     
    5457#  include "vectopt_loop_substitute.h90" 
    5558   !!---------------------------------------------------------------------- 
    56    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     59   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    5760   !! $Id$ 
    5861   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    6871      !!      the temporal evolution of some key variables 
    6972      !!------------------------------------------------------------------- 
    70       INTEGER, INTENT(in) ::   kt   ! number of iteration 
     73      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    7174      !! 
    7275      INTEGER  ::   jv,ji, jj   ! dummy loop indices 
    7376      INTEGER  ::   nv          ! indice of variable  
    74       REAL(wp) ::   zarea    , zldarea     ! sea-ice and leads area 
    75       REAL(wp) ::   zextent15, zextent85   ! sea-ice extent (15% and 85%) 
    76       REAL(wp) ::   zicevol  , zsnwvol     ! sea-ice and snow volume volume 
    77       REAL(wp) ::   zicespd                     ! sea-ice velocity 
     77      REAL(wp) ::   zarea    , zldarea  ,    &  ! sea-ice and leads area 
     78         &          zextent15, zextent85,    &  ! sea-ice extent (15% and 85%) 
     79         &          zicevol  , zsnwvol  ,    &  ! sea-ice and snow volume volume 
     80         &          zicespd                     ! sea-ice velocity 
    7881      REAL(wp), DIMENSION(jpinfmx) ::   vinfor  ! temporary working space  
    7982      !!------------------------------------------------------------------- 
     
    117120      vinfor(nv+13) = SQRT( vinfor(nv+13) / MAX( vinfor(nv+9) , epsi06 ) ) 
    118121 
     122 
    119123      ! variables in southern Hemis 
    120124       nv = nv + 1 
     
    173177       INTEGER  ::   nv            ! indice of variable  
    174178       REAL(wp) ::   zxx0, zxx1    ! temporary scalars 
    175        !! 
     179 
    176180       NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy 
    177181       !!------------------------------------------------------------------- 
    178182 
    179        REWIND ( numnam_ice )                 ! Read Namelist namicedia 
     183       ! Read Namelist namicedia 
     184       REWIND ( numnam_ice ) 
    180185       READ   ( numnam_ice  , namicedia ) 
    181186        
    182        IF(lwp) THEN                          ! control print 
     187       IF(lwp) THEN 
    183188          WRITE(numout,*) 
    184189          WRITE(numout,*) 'lim_dia_init_2 : ice parameters for ice diagnostics ' 
     
    269274   !!   Default option :                           NO LIM 2.0 sea-ice model 
    270275   !!---------------------------------------------------------------------- 
     276CONTAINS 
     277   SUBROUTINE lim_dia_2         ! Empty routine 
     278   END SUBROUTINE lim_dia_2 
    271279#endif 
    272280 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdmp_2.F90

    r1855 r1857  
    44   !!  Ice model : restoring Ice thickness and Fraction leads 
    55   !!====================================================================== 
    6    !! History :   2.0  ! 2004-04 (S. Theetten) Original code 
     6   !! History :   2.0  !  04-04 (S. Theetten) Original code 
    77   !!---------------------------------------------------------------------- 
    88#if defined key_lim2   &&   defined key_tradmp 
     
    1111   !!   'key_tradmp'                                                Damping 
    1212   !!---------------------------------------------------------------------- 
    13    !!   lim_dmp_2       : ice model damping 
    1413   !!---------------------------------------------------------------------- 
    15    USE oce              ! ocean variables 
    16    USE dom_oce          ! ocean domain 
    17    USE phycst           ! physical constants 
    18    USE ice_2            ! LIM-2 variables 
    19    USE tradmp           ! traceur damping 
    20    USE in_out_manager   ! I/O manager 
    21    USE iom              ! 
     14   !!   lim_dmp_2      : ice model damping 
     15   !!---------------------------------------------------------------------- 
     16   USE in_out_manager  ! I/O manager 
     17   USE phycst          ! physical constants 
     18   USE ice_2 
     19   USE tradmp 
     20   USE dom_oce 
     21   USE oce 
     22   USE iom 
    2223    
    2324   IMPLICIT NONE 
     
    2627   PUBLIC   lim_dmp_2     ! called by ice_step_2 
    2728    
    28    INTEGER                        ::   nice1, nice2   ! first and second record used 
    29    INTEGER                        ::   inumice_dmp    ! logical unit for ice variables (damping) 
    30    REAL(wp), DIMENSION(jpi,jpj)   ::   hicif_dta      ! ice thickness at a given time 
    31    REAL(wp), DIMENSION(jpi,jpj)   ::   frld_dta       ! fraction lead at a given time 
    32    REAL(wp), DIMENSION(jpi,jpj,2) ::   hicif_data     ! ice thickness data at two consecutive times 
    33    REAL(wp), DIMENSION(jpi,jpj,2) ::   frld_data      ! fraction lead data at two consecutive times 
     29   INTEGER                        ::   nice1, nice2,  &  ! first and second record used 
     30      &                                inumice_dmp       ! logical unit for ice variables (damping) 
     31   REAL(wp), DIMENSION(jpi,jpj)   ::   hicif_dta  ,   &  ! ice thickness at a given time 
     32      &                                frld_dta          ! fraction lead at a given time 
     33   REAL(wp), DIMENSION(jpi,jpj,2) ::   hicif_data ,   &  ! ice thickness data at two consecutive times 
     34      &                                frld_data         ! fraction lead data at two consecutive times 
    3435 
    3536   !! * Substitution 
    3637#  include "vectopt_loop_substitute.h90" 
    3738   !!---------------------------------------------------------------------- 
    38    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     39   !!   LIM 2.0 , UCL-LOCEAN-IPSL  (2006) 
    3940   !! $Id$ 
    4041   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    4344CONTAINS 
    4445 
    45    SUBROUTINE lim_dmp_2( kt ) 
     46   SUBROUTINE lim_dmp_2(kt) 
    4647      !!------------------------------------------------------------------- 
    4748      !!                   ***  ROUTINE lim_dmp_2 *** 
     
    5253      !! ** method  : the key_tradmp must be used to compute resto(:,:) coef. 
    5354      !!--------------------------------------------------------------------- 
    54       INTEGER, INTENT(in) ::   kt   ! ocean time-step 
    55       !! 
    56       INTEGER ::   ji, jj   ! dummy loop indices 
     55      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
     56      ! 
     57      INTEGER             ::   ji, jj         ! dummy loop indices 
    5758      !!--------------------------------------------------------------------- 
    5859      ! 
    5960      CALL dta_lim_2( kt ) 
    60       ! 
     61 
    6162      DO jj = 2, jpjm1 
    6263         DO ji = fs_2, fs_jpim1   ! vector opt. 
    6364            hicif(ji,jj) = hicif(ji,jj) - rdt_ice * resto(ji,jj,1) * ( hicif(ji,jj) - hicif_dta(ji,jj) ) 
    64             frld(ji,jj)  = frld (ji,jj) - rdt_ice * resto(ji,jj,1) * ( frld (ji,jj) - frld_dta (ji,jj) )   
     65            frld(ji,jj)  = frld (ji,jj) - rdt_ice * resto(ji,jj,1) * ( frld(ji,jj) - frld_dta (ji,jj) )   
    6566         END DO 
    6667      END DO 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdyn_2.F90

    r1855 r1857  
    1717   !!---------------------------------------------------------------------- 
    1818   USE dom_oce        ! ocean space and time domain 
    19    USE sbc_oce        ! surface boundary condition : ocean 
    20    USE phycst         ! physical constants 
    21    USE ice_2          ! LIM-2 sea-ice variables 
    22    USE dom_ice_2      ! LIM-2 domain 
    23    USE limistate_2    ! LIM-2 initial state 
    24    USE limrhg_2       ! LIM-2 rheology 
     19   USE sbc_oce        ! 
     20   USE phycst         ! 
     21   USE ice_2          ! 
     22   USE dom_ice_2      ! 
     23   USE limistate_2    ! 
     24   USE limrhg_2       ! ice rheology 
    2525 
    2626   USE lbclnk         ! 
     
    3434   PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim 
    3535 
     36   !! * Module variables 
    3637   REAL(wp)  ::  rone    = 1.e0   ! constant value 
    3738 
    3839#  include "vectopt_loop_substitute.h90" 
    3940   !!---------------------------------------------------------------------- 
    40    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     41   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
    4142   !! $Id$ 
    4243   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    5859      !!              - treatment of the case if no ice dynamic 
    5960      !!--------------------------------------------------------------------- 
    60       INTEGER, INTENT(in) ::   kt   ! number of iteration 
    61       !! 
    62       INTEGER  ::   ji, jj        ! dummy loop indices 
    63       INTEGER  ::   i_j1, i_jpj   ! Starting/ending j-indices for rheology 
    64       REAL(wp) ::   zcoef         ! temporary scalar 
     61      INTEGER, INTENT(in) ::   kt     ! number of iteration 
     62      !! 
     63      INTEGER  ::   ji, jj             ! dummy loop indices 
     64      INTEGER  ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology 
     65      REAL(wp) ::   zcoef              ! temporary scalar 
    6566      REAL(wp), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice 
    6667      REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask 
     
    143144         ENDIF 
    144145 
    145          IF(ln_ctl)   CALL prt_ctl( tab2d_1=u_ice, clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :' ) 
     146         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 
    146147          
    147148         ! computation of friction velocity 
     
    178179      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point 
    179180      ! 
    180       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ust2s, clinfo1=' lim_dyn  : ust2s :' ) 
     181      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    181182 
    182183   END SUBROUTINE lim_dyn_2 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limhdf_2.F90

    r1855 r1857  
    44   !! LIM 2.0 ice model : horizontal diffusion of sea-ice quantities 
    55   !!====================================================================== 
    6    !! History :  LIM  ! 2000-01 (UCL)  Original code 
    7    !!            1.0  ! 2001-05 (G. Madec, R. Hordoir) opa norm 
    8    !!            2.0  ! 2002-08 (C. Ethe)  F90, free form 
    9    !!---------------------------------------------------------------------- 
    106#if defined key_lim2 
    117   !!---------------------------------------------------------------------- 
     
    1410   !!   lim_hdf_2  : diffusion trend on sea-ice variable 
    1511   !!---------------------------------------------------------------------- 
    16    USE dom_oce          ! ocean domain 
    17    USE ice_2            ! LIM-2 variables 
    18    USE prtctl           ! Print control 
    19    USE lbclnk           ! 
    20    USE lib_mpp          ! 
    21    USE in_out_manager   ! I/O manager 
     12   !! * Modules used 
     13   USE dom_oce 
     14   USE in_out_manager 
     15   USE ice_2 
     16   USE lbclnk 
     17   USE lib_mpp 
     18   USE prtctl          ! Print control 
    2219 
    2320   IMPLICIT NONE 
    2421   PRIVATE 
    2522 
    26    PUBLIC   lim_hdf_2   ! called by lim_tra_2 
     23   !! * Routine accessibility 
     24   PUBLIC lim_hdf_2    ! called by lim_tra_2 
    2725 
    28    LOGICAL  ::   linit = .TRUE.              ! indictor of initialisation 
     26   !! * Module variables 
     27   LOGICAL  ::   linit = .TRUE.              ! ??? 
    2928   REAL(wp) ::   epsi04 = 1e-04              ! constant 
    30    REAL(wp), DIMENSION(jpi,jpj) ::   zfact   ! metric term 
     29   REAL(wp), DIMENSION(jpi,jpj) ::   zfact   ! ??? 
    3130 
    3231   !! * Substitution  
    3332#  include "vectopt_loop_substitute.h90" 
    3433   !!---------------------------------------------------------------------- 
    35    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     34   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    3635   !! $Id$ 
    37    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     36   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3837   !!---------------------------------------------------------------------- 
    3938 
     
    4443      !!                  ***  ROUTINE lim_hdf_2  *** 
    4544      !! 
    46       !! ** purpose :   Compute and add the diffusive trend on sea-ice variables 
     45      !! ** purpose :   Compute and add the diffusive trend on sea-ice 
     46      !!      variables 
    4747      !! 
    4848      !! ** method  :   Second order diffusive operator evaluated using a 
    49       !!              Cranck-Nicholson time Scheme. 
     49      !!      Cranck-Nicholson time Scheme. 
    5050      !! 
    51       !! ** Action  :   update ptab with the diffusive contribution 
     51      !! ** Action  :    update ptab with the diffusive contribution 
     52      !! 
     53      !! History : 
     54      !!        !  00-01 (LIM) Original code 
     55      !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
     56      !!        !  02-08 (C. Ethe)  F90, free form 
    5257      !!------------------------------------------------------------------- 
    53       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   ptab   ! Field on which the diffusion is applied   
    54       !! 
     58      ! * Arguments 
     59      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   & 
     60         ptab                 ! Field on which the diffusion is applied   
     61      REAL(wp), DIMENSION(jpi,jpj) ::   & 
     62         ptab0                ! ??? 
     63 
     64      ! * Local variables 
    5565      INTEGER ::  ji, jj      ! dummy loop indices 
    56       INTEGER ::  its, iter   ! temporary integers 
     66      INTEGER ::  & 
     67         its, iter            ! temporary integers 
    5768      CHARACTER (len=55) :: charout 
    58       REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! temporary scalars 
    59       REAL(wp), DIMENSION(jpi,jpj) ::   zrlx, zflu, zflv, zdiv0, zdiv  ! temporary workspaces 
    60       REAL(wp), DIMENSION(jpi,jpj) ::   ztab0                ! ??? 
     69      REAL(wp) ::  & 
     70         zalfa, zrlxint, zconv, zeps   ! temporary scalars 
     71      REAL(wp), DIMENSION(jpi,jpj) ::  &  
     72         zrlx, zflu, zflv, &  ! temporary workspaces 
     73         zdiv0, zdiv          !    "           " 
    6174      !!------------------------------------------------------------------- 
    6275 
     
    6982 
    7083      ! Arrays initialization 
    71       ztab0(:, : ) = ptab(:,:) 
     84      ptab0 (:, : ) = ptab(:,:) 
     85!bug  zflu (:,jpj) = 0.e0 
     86!bug  zflv (:,jpj) = 0.e0 
    7287      zdiv0(:, 1 ) = 0.e0 
    7388      zdiv0(:,jpj) = 0.e0 
     
    8398         DO jj = 2, jpjm1   
    8499            DO ji = fs_2 , fs_jpim1   ! vector opt. 
    85                zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
     100               zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj  ) + e1v(ji,jj) + e1v(ji,jj-1) ) & 
     101                  &          / ( e1t(ji,jj) * e2t(ji,jj) ) 
    86102            END DO 
    87103         END DO 
     
    94110      iter  = 0 
    95111 
    96       !                                                   !======================! 
    97       DO WHILE ( ( zconv > zeps ) .AND. (iter <= its) )   !  Sub-time step loop  ! 
    98          !                                                !======================! 
    99          iter = iter + 1               ! incrementation of the sub-time step number 
     112      !                                                   !=================== 
     113      DO WHILE ( ( zconv > zeps ) .AND. (iter <= its) )   ! Sub-time step loop 
     114         !                                                !=================== 
     115         ! incrementation of the sub-time step number 
     116         iter = iter + 1 
    100117 
    101          DO jj = 1, jpjm1              ! diffusive fluxes in U- and V- direction 
     118         ! diffusive fluxes in U- and V- direction 
     119         DO jj = 1, jpjm1 
    102120            DO ji = 1 , fs_jpim1   ! vector opt. 
    103121               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
     
    105123            END DO 
    106124         END DO 
    107          DO jj= 2, jpjm1               ! diffusive trend : divergence of the fluxes 
     125 
     126         ! diffusive trend : divergence of the fluxes 
     127         DO jj= 2, jpjm1 
    108128            DO ji = fs_2 , fs_jpim1   ! vector opt.  
    109129               zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   & 
     
    111131            END DO 
    112132         END DO 
    113          !                             ! save the first evaluation of the diffusive trend in zdiv0 
     133 
     134         ! save the first evaluation of the diffusive trend in zdiv0 
    114135         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        
    115136 
    116          DO jj = 2, jpjm1              ! XXXX iterative evaluation????? 
     137         ! XXXX iterative evaluation????? 
     138         DO jj = 2, jpjm1 
    117139            DO ji = fs_2 , fs_jpim1   ! vector opt. 
    118                zrlxint = (   ztab0(ji,jj)    & 
     140               zrlxint = (   ptab0(ji,jj)    & 
    119141                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + zfact(ji,jj) * ptab(ji,jj) )   & 
    120142                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             &  
     
    123145            END DO 
    124146         END DO 
     147 
     148         ! lateral boundary condition on ptab 
    125149         CALL lbc_lnk( zrlx, 'T', 1. ) 
    126150 
    127          zconv = 0.e0                  ! convergence test 
     151         ! convergence test 
     152         zconv = 0.e0 
    128153         DO jj = 2, jpjm1 
    129154            DO ji = 2, jpim1 
     
    132157         END DO 
    133158         IF( lk_mpp )   CALL mpp_max( zconv )   ! max over the global domain 
    134          ! 
    135          ptab(:,:) = zrlx(:,:)         ! update value 
    136          !                                         !=============================! 
    137       END DO                                       !  end of sub-time step loop  ! 
    138       !                                            !=============================! 
     159 
     160         ptab(:,:) = zrlx(:,:) 
     161 
     162         !                                         !========================== 
     163      END DO                                       ! end of sub-time step loop 
     164      !                                            !========================== 
    139165 
    140166      IF(ln_ctl)   THEN 
    141          zrlx(:,:) = ptab(:,:) - ztab0(:,:) 
     167         zrlx(:,:) = ptab(:,:) - ptab0(:,:) 
    142168         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter 
    143          CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout ) 
     169         CALL prt_ctl(tab2d_1=zrlx, clinfo1=charout) 
    144170      ENDIF 
    145       ! 
     171 
    146172   END SUBROUTINE lim_hdf_2 
    147173 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limistate_2.F90

    r1855 r1857  
    44   !!              Initialisation of diagnostics ice variables 
    55   !!====================================================================== 
    6    !! History :   1.0  !  2001-04  (C. Ethe, G. Madec)  Original code 
    7    !!             2.0  !  2003-08  (G. Madec)  add lim_istate_init 
    8    !!              -   !  2004-04  (S. Theetten) initialization from a file 
    9    !!              -   !  2006-07  (S. Masson)  IOM to read the restart 
    10    !!              -   !  2007-10  (G. Madec)  surface module 
     6   !! History :   1.0  !  01-04  (C. Ethe, G. Madec)  Original code 
     7   !!             2.0  !  03-08  (G. Madec)  add lim_istate_init 
     8   !!                  !  04-04  (S. Theetten) initialization from a file 
     9   !!                  !  06-07  (S. Masson)  IOM to read the restart 
     10   !!                  !  07-10  (G. Madec)  surface module 
    1111   !!-------------------------------------------------------------------- 
    1212#if defined key_lim2 
     
    1414   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1515   !!---------------------------------------------------------------------- 
     16   !!---------------------------------------------------------------------- 
    1617   !!   lim_istate_2      :  Initialisation of diagnostics ice variables 
    1718   !!   lim_istate_init_2 :  initialization of ice state and namelist read 
    1819   !!---------------------------------------------------------------------- 
    19    USE oce              ! ocean variables 
    20    USE ice_2            ! LIM-2 variables 
    21    USE par_ice_2        ! LIM-2 ice parameters 
    22    USE dom_ice_2        ! LIM-2 domain 
    23    USE phycst           ! physical constants 
    24    USE eosbn2           ! equation of state 
    25    USE lbclnk           !  
    26    USE iom              ! 
    27    USE in_out_manager   ! 
     20   USE phycst 
     21   USE par_ice_2       ! ice parameters 
     22   USE dom_ice_2 
     23   USE eosbn2          ! equation of state 
     24   USE lbclnk 
     25   USE oce 
     26   USE ice_2 
     27   USE iom 
     28   USE in_out_manager 
    2829 
    2930   IMPLICIT NONE 
    3031   PRIVATE 
    3132 
    32    PUBLIC   lim_istate_2   ! routine called by lim_init_2.F90 
    33  
    34    !                                 !!! ** init namelist (namiceini) ** 
    35    LOGICAL  ::   ln_limini = .FALSE.  ! Ice initialization state 
     33   PUBLIC lim_istate_2      ! routine called by lim_init_2.F90 
     34 
     35!!! ** init namelist (namiceini) ** 
     36   LOGICAL  ::   ln_limini = .FALSE.  !: Ice initialization state 
    3637   REAL(wp) ::   ttest     = 2.0      ! threshold water temperature for initial sea ice 
    3738   REAL(wp) ::   hninn     = 0.5      ! initial snow thickness in the north 
     
    4445   REAL(wp) ::   zero      = 0.e0     ! constant value = 0 
    4546   REAL(wp) ::   zone      = 1.e0     ! constant value = 1 
    46     
    47    !!---------------------------------------------------------------------- 
    48    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     47   !!---------------------------------------------------------------------- 
     48   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
    4949   !! $Id$ 
    5050   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    6262      !!                or from arbitrary sea-ice conditions 
    6363      !!-------------------------------------------------------------------- 
    64       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    65       REAL(wp) ::   zidto        ! temporary scalar 
     64      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     65      REAL(wp) ::   zidto                     ! temporary scalar 
    6666      !-------------------------------------------------------------------- 
    6767  
     
    6969 
    7070      IF( .NOT. ln_limini ) THEN   
    71          ! 
     71          
    7272         tfu(:,:) = tfreez( sn(:,:,1) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
    73          ! 
     73 
    7474         DO jj = 1, jpj 
    7575            DO ji = 1, jpi 
     
    126126 
    127127      !-- lateral boundary conditions 
    128       CALL lbc_lnk( hicif, 'T', 1. )   ;   CALL lbc_lnk( frld , 'T', 1. ) 
     128      CALL lbc_lnk( hicif, 'T', 1. ) 
     129      CALL lbc_lnk( frld , 'T', 1. ) 
    129130 
    130131      ! C A U T I O N  frld = 1 over land and lbc_lnk put zero along  
     
    139140      CALL lbc_lnk( fsbbq  , 'T', 1. ) 
    140141      CALL lbc_lnk( qstoif , 'T', 1. ) 
    141       ! 
     142 
    142143   END SUBROUTINE lim_istate_2 
    143144 
     
    150151      !! 
    151152      !! ** Method  :   Read the namiceini namelist and check the parameter  
    152       !!              values called at the first timestep (nit000) 
     153      !!       values called at the first timestep (nit000) 
    153154      !! 
    154155      !! ** input   :   Namelist namiceini 
    155156      !!------------------------------------------------------------------- 
    156       INTEGER ::   ji,jj      ! dummy loop indices 
    157       INTEGER ::   inum_ice   ! temporary integer 
    158       !! 
     157      INTEGER :: inum_ice 
     158      INTEGER :: ji,jj 
     159 
    159160      NAMELIST/namiceini/ ln_limini, ttest, hninn, hginn, alinn, & 
    160          &                                  hnins, hgins, alins 
     161         &                hnins, hgins, alins 
    161162      !!------------------------------------------------------------------- 
    162163      ! 
     
    164165      READ   ( numnam_ice , namiceini ) 
    165166      ! 
    166       IF(lwp) THEN                        ! control print 
     167      IF(lwp) THEN 
    167168         WRITE(numout,*) 
    168169         WRITE(numout,*) 'lim_istate_init_2 : ice parameters inititialisation ' 
     
    178179      ENDIF 
    179180 
    180       IF( ln_limini ) THEN                ! Ice initialization using input file 
     181      IF( ln_limini ) THEN                      ! Ice initialization using input file 
    181182         ! 
    182183         CALL iom_open( 'Ice_initialization.nc', inum_ice ) 
     
    185186            IF(lwp) WRITE(numout,*) 
    186187            IF(lwp) WRITE(numout,*) '                  ice state initialization with : Ice_initialization.nc' 
    187             ! 
     188             
    188189            CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif )       
    189190            CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif )       
     
    191192            CALL iom_get( inum_ice, jpdom_data, 'ts'   , sist  ) 
    192193            CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif(1:nlci,1:nlcj,:),   & 
    193                        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) ) 
     194                 &        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) ) 
    194195            ! put some values in the extra-halo... 
    195196            DO jj = nlcj+1, jpj   ;   tbif(1:nlci,jj,:) = tbif(1:nlci,nlej,:)   ;   END DO 
    196197            DO ji = nlci+1, jpi   ;   tbif(ji    ,: ,:) = tbif(nlei  ,:   ,:)   ;   END DO 
    197             ! 
     198 
    198199            CALL iom_close( inum_ice) 
    199200            ! 
     
    207208   !!   Default option :         Empty module      NO LIM 2.0 sea-ice model 
    208209   !!---------------------------------------------------------------------- 
     210CONTAINS 
     211   SUBROUTINE lim_istate_2        ! Empty routine 
     212   END SUBROUTINE lim_istate_2 
    209213#endif 
    210214 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limmsh_2.F90

    r1855 r1857  
    44   !! LIM 2.0 ice model :   definition of the ice mesh parameters 
    55   !!====================================================================== 
    6    !! ** History :  LIM  ! 2001-04 (UCL)  original code 
    7    !!               2.0  ! 2002-08 (C. Ethe, G. Madec) F90, module 
    8    !!---------------------------------------------------------------------- 
    96#if defined key_lim2 
    107   !!---------------------------------------------------------------------- 
     
    1310   !!   lim_msh_2   : definition of the ice mesh 
    1411   !!---------------------------------------------------------------------- 
     12   !! * Modules used 
    1513   USE phycst 
    1614   USE dom_oce 
     
    2220   PRIVATE 
    2321 
    24    PUBLIC   lim_msh_2   ! routine called by ice_ini_2.F90 
    25  
    26    !!---------------------------------------------------------------------- 
    27    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     22   !! * Accessibility 
     23   PUBLIC lim_msh_2      ! routine called by ice_ini_2.F90 
     24 
     25   !!---------------------------------------------------------------------- 
     26   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    2827   !! $Id$ 
    29    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     28   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3029   !!---------------------------------------------------------------------- 
    3130 
     
    4342      !!              - Initialization of the ice masks (tmsk, umsk) 
    4443      !!  
    45       !! Reference  : Deleersnijder et al. Ocean Modelling 100, 7-10  
     44      !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10  
     45      !! 
     46      !! ** History : 
     47      !!         original    : 01-04 (LIM) 
     48      !!         addition    : 02-08 (C. Ethe, G. Madec) 
    4649      !!---------------------------------------------------------------------  
    47       INTEGER  ::   ji, jj   ! dummy loop indices 
    48       REAL(wp) ::   zh1p, zd1d2p, zusden    ! temporary scalars 
    49       REAL(wp) ::   zh2p, zd2d1p, zusden2   !    -         - 
    50       REAL(wp), DIMENSION(jpi,jpj) ::   zd2d1 , zd1d2   ! 2D workspace 
     50      !! * Local variables 
     51      INTEGER :: ji, jj      ! dummy loop indices 
     52 
     53      REAL(wp), DIMENSION(jpi,jpj) ::  & 
     54         zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction 
     55         !                   ! (resp. y direction) (defined at the center) 
     56      REAL(wp) ::         & 
     57         zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid 
     58         zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid 
     59         zusden, zusden2     ! temporary scalars 
    5160      !!--------------------------------------------------------------------- 
    5261 
     
    103112      !------------------- 
    104113!!ibug ??? 
    105       akappa(:,:,:,:)     = 0.e0 
    106       wght  (:,:,:,:)    = 0.e0 
     114      akappa(:,:,:,:) = 0.e0 
     115      wght(:,:,:,:) = 0.e0 
    107116      alambd(:,:,:,:,:,:) = 0.e0 
    108       tmu   (:,:)        = 0.e0 
     117      tmu(:,:) = 0.e0 
    109118!!i 
    110119       
     
    229238         END DO 
    230239      END DO 
    231       CALL lbc_lnk( tmu(:,:), 'I', 1. )      !--lateral boundary conditions     
    232        
    233       area(:,:) = e1t(:,:) * e2t(:,:)        ! unmasked and masked area of T-grid cell 
     240       
     241      !--lateral boundary conditions     
     242      CALL lbc_lnk( tmu(:,:), 'I', 1. ) 
     243       
     244      ! unmasked and masked area of T-grid cell 
     245      area(:,:) = e1t(:,:) * e2t(:,:) 
    234246       
    235247   END SUBROUTINE lim_msh_2 
     
    239251   !!   Default option            Dummy Module         NO LIM sea-ice model 
    240252   !!---------------------------------------------------------------------- 
     253CONTAINS 
     254   SUBROUTINE lim_msh_2           ! Dummy routine 
     255   END SUBROUTINE lim_msh_2 
    241256#endif 
    242257 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limrhg_2.F90

    r1855 r1857  
    44   !!   Ice rheology :  performs sea ice rheology 
    55   !!====================================================================== 
    6    !! History :  0.0  !  1993-12  (M.A. Morales Maqueda.)  Original code 
    7    !!            1.0  !  1994-12  (H. Goosse)  
    8    !!            2.0  !  1903-12  (C. Ethe, G. Madec)  F90, mpp 
    9    !!             -   !  2006-08  (G. Madec)  surface module, ice-stress at I-point 
    10    !!             -   !  2009-09  (G. Madec)  Huge verctor optimisation 
     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_lim2 
     
    1414   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    1515   !!---------------------------------------------------------------------- 
    16    !!   lim_rhg_2     : computes ice velocities 
     16   !!---------------------------------------------------------------------- 
     17   !!   lim_rhg_2   : computes ice velocities 
    1718   !!---------------------------------------------------------------------- 
    1819   USE par_oce        ! ocean parameter 
     
    3334   PUBLIC   lim_rhg_2 ! routine called by lim_dyn 
    3435 
    35    REAL(wp) ::   rzero = 0.e0   ! constant value: zero 
    36    REAL(wp) ::   rone  = 1.e0   !            and  one 
     36   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero 
     37   REAL(wp) ::   rone    = 1.e0   !            and  one 
    3738 
    3839   !! * Substitutions 
    3940#  include "vectopt_loop_substitute.h90" 
    4041   !!---------------------------------------------------------------------- 
    41    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     42   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
    4243   !! $Id$ 
    4344   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    5556      !!  viscous-plastic law including shear strength and a bulk rheology. 
    5657      !! 
    57       !! ** Action  : - compute u_ice, v_ice the sea-ice velocity defined at I-point 
     58      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity defined 
     59      !!              at I-point 
    5860      !!------------------------------------------------------------------- 
    5961      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation 
     
    580582   !!   Default option          Dummy module       NO 2.0 LIM sea-ice model 
    581583   !!---------------------------------------------------------------------- 
     584CONTAINS 
     585   SUBROUTINE lim_rhg_2( k1 , k2 )       ! Dummy routine 
     586      WRITE(*,*) 'lim_rhg_2: You should not have seen this print! error?', k1, k2 
     587   END SUBROUTINE lim_rhg_2 
    582588#endif 
    583589 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limrst_2.F90

    r1855 r1857  
    22   !!====================================================================== 
    33   !!                     ***  MODULE  limrst_2  *** 
    4    !! LIM-2 ice model : open/read/write the restart file 
     4   !! Ice restart :  write the ice restart file 
    55   !!====================================================================== 
    6    !! History :  2.0  ! 2001-04  (C. Ethe, G. Madec)  Original code 
    7    !!                 ! 2006-07  (S. Masson)  use IOM for restart read/write 
     6   !! History :  2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
     7   !!                 !  06-07  (S. Masson)  use IOM for restart read/write 
    88   !!---------------------------------------------------------------------- 
    99#if defined key_lim2 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
     12   !!---------------------------------------------------------------------- 
    1213   !!---------------------------------------------------------------------- 
    1314   !!   lim_rst_opn_2   : open ice restart file 
     
    1920   USE sbc_oce 
    2021   USE sbc_ice 
     22 
    2123   USE in_out_manager 
    2224   USE iom 
     
    3335 
    3436   !!---------------------------------------------------------------------- 
    35    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     37   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
    3638   !! $Id$ 
    3739   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    8284   END SUBROUTINE lim_rst_opn_2 
    8385 
    84  
    8586   SUBROUTINE lim_rst_write_2( kt ) 
    8687      !!---------------------------------------------------------------------- 
     
    9091      !!---------------------------------------------------------------------- 
    9192      INTEGER, INTENT(in) ::   kt   ! number of iteration 
    92       !! 
     93      ! 
    9394      INTEGER ::   iter   ! kt + nn_fsbc -1 
    9495      !!---------------------------------------------------------------------- 
     
    180181      ENDIF 
    181182 
    182       IF( jprstlib == jprstdimg ) THEN 
     183      IF ( jprstlib == jprstdimg ) THEN 
    183184        ! eventually read netcdf file (monobloc)  for restarting on different number of processors 
    184185        ! if {cn_icerst_in}.nc exists, then set jlibalt to jpnf90 
     
    187188      ENDIF 
    188189 
    189       CALL iom_open( cn_icerst_in, numrir, kiolib = jlibalt ) 
     190      CALL iom_open ( cn_icerst_in, numrir, kiolib = jlibalt ) 
    190191 
    191192      CALL iom_get( numrir, 'kt_ice' , ziter ) 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90

    r1855 r1857  
    22   !!====================================================================== 
    33   !!                       ***  MODULE limsbc_2   *** 
    4    !!   LIM 2 ice model :   heat, salt, mass and momentum fluxes at the sea ice/ocean interface 
     4   !!           computation of the flux at the sea ice/ocean interface 
    55   !!====================================================================== 
    6    !! History :  LIM  ! 2000-01 (H. Goosse) Original code 
    7    !!            2.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90 
    8    !!             -   ! 2006-07 (G. Madec) surface module 
     6   !! History : 00-01 (H. Goosse) Original code 
     7   !!           02-07 (C. Ethe, G. Madec) re-writing F90 
     8   !!           06-07 (G. Madec) surface module 
    99   !!---------------------------------------------------------------------- 
    1010#if defined key_lim2 
    1111   !!---------------------------------------------------------------------- 
    1212   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
     13   !!---------------------------------------------------------------------- 
    1314   !!---------------------------------------------------------------------- 
    1415   !!   lim_sbc_2  : flux at the ice / ocean interface 
     
    1920   USE sbc_oce          ! surface boundary condition 
    2021   USE phycst           ! physical constants 
    21    USE ice_2            ! LIM 2 sea-ice variables 
     22   USE ice_2            ! LIM sea-ice variables 
    2223 
    2324   USE lbclnk           ! ocean lateral boundary condition 
     
    3233   PRIVATE 
    3334 
    34    PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2 
     35   PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2 
    3536 
    3637   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values 
     
    4344#  include "vectopt_loop_substitute.h90" 
    4445   !!---------------------------------------------------------------------- 
    45    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     46   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
    4647   !! $Id$ 
    4748   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    108109         sice_r(:,:) = sice 
    109110         ! 
    110          IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN             !  ORCA_R2 configuration 
     111         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN 
     112            !                                        ! ======================= 
     113            !                                        !  ORCA_R2 configuration 
     114            !                                        ! ======================= 
    111115            ii0 = 145   ;   ii1 = 180        ! Baltic Sea 
    112116            ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 
     
    190194      END DO 
    191195 
    192       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:)                         )       
    193       CALL iom_put( 'qns_io_cea'  , qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )       
    194       CALL iom_put( 'qsr_io_cea'  , fstric(:,:) * (1. - pfrld(:,:))      ) 
     196      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )       
     197      CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )       
     198      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) ) 
    195199 
    196200      !------------------------------------------! 
     
    298302      !-----------------------------------------------! 
    299303 
    300       IF( lk_cpl ) THEN            
     304      IF ( lk_cpl ) THEN            
     305         ! Ice surface temperature  
    301306         tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature        
    302          !                                  ! snow/ice and ocean albedos 
     307         ! Computation of snow/ice and ocean albedo 
    303308         CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb ) 
    304309         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys) 
     
    307312 
    308313      IF(ln_ctl) THEN 
    309          CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns  , clinfo2=' qns     : ' ) 
    310          CALL prt_ctl( tab2d_1=emp , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps , clinfo2=' emps    : ' ) 
    311          CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   & 
    312             &          tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask ) 
    313          CALL prt_ctl( tab2d_1=fr_i, clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ' ) 
     314         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ') 
     315         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ') 
     316         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   & 
     317            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask ) 
     318         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ') 
    314319      ENDIF  
    315       ! 
     320    
     321    END SUBROUTINE lim_sbc_2 
     322 
     323#else 
     324   !!---------------------------------------------------------------------- 
     325   !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model 
     326   !!---------------------------------------------------------------------- 
     327CONTAINS 
     328   SUBROUTINE lim_sbc_2         ! Dummy routine 
    316329   END SUBROUTINE lim_sbc_2 
    317  
    318 #else 
    319    !!---------------------------------------------------------------------- 
    320    !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model 
    321    !!---------------------------------------------------------------------- 
    322330#endif  
    323331 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtab_2.F90

    r1855 r1857  
    22   !!====================================================================== 
    33   !!                       ***  MODULE limtab_2   *** 
    4    !!   LIM 2 ice model : transform 1D (2D) array to a 2D (1D) array 
     4   !!              transform 1D (2D) array to a 2D (1D) table 
    55   !!====================================================================== 
    66#if defined key_lim2 
    77   !!---------------------------------------------------------------------- 
    8    !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
     8   !!   tab_2d_1d  : 2-D to 1-D 
     9   !!   tab_1d_2d  : 1-D to 2-D 
    910   !!---------------------------------------------------------------------- 
    10    !!   tab_2d_1d_2  : 2-D to 1-D 
    11    !!   tab_1d_2d_2  : 1-D to 2-D 
    12    !!---------------------------------------------------------------------- 
     11   !! * Modules used 
    1312   USE par_kind 
    1413 
     
    1615   PRIVATE 
    1716 
    18    PUBLIC   tab_2d_1d_2   ! called by lim_thd_2 
    19    PUBLIC   tab_1d_2d_2   ! called by lim_thd_2 
     17   !! * Routine accessibility 
     18   PUBLIC tab_2d_1d_2  ! called by lim_ther 
     19   PUBLIC tab_1d_2d_2  ! called by lim_ther 
    2020 
    2121   !!---------------------------------------------------------------------- 
    22    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     22   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    2323   !! $Id$ 
    2424   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    2626CONTAINS 
    2727 
    28    SUBROUTINE tab_2d_1d_2( kdim1d, ptab1d, ptab2d, kdim2d_i, kdim2d_j, ptab_ind ) 
    29       !!------------------------------------------------------------------- 
    30       INTEGER , INTENT(in   )                                ::   kdim1d 
    31       INTEGER , INTENT(in   )                                ::   kdim2d_i, kdim2d_j 
    32       REAL(wp), INTENT(in   ), DIMENSION(kdim2d_i, kdim2d_j) ::   ptab2d 
    33       INTEGER , INTENT(in   ), DIMENSION(kdim1d)             ::   ptab_ind 
    34       REAL(wp), INTENT(  out), DIMENSION(kdim1d)             ::   ptab1d 
    35       !! 
    36       INTEGER ::   jn , jid, jjd   ! dummy loop indices 
    37       !!------------------------------------------------------------------- 
    38       !  
    39       DO jn = 1, kdim1d 
    40          jid = MOD( ptab_ind(jn) - 1, kdim2d_i ) + 1 
    41          jjd = ( ptab_ind(jn) - 1 ) / kdim2d_i + 1 
    42          ptab1d(jn) = ptab2d(jid,jjd) 
     28   SUBROUTINE tab_2d_1d_2 ( ndim1d, tab1d, tab2d, ndim2d_x, ndim2d_y, tab_ind ) 
     29 
     30      INTEGER, INTENT(in) :: & 
     31         ndim1d, ndim2d_x, ndim2d_y 
     32 
     33      REAL(wp), DIMENSION (ndim2d_x, ndim2d_y), INTENT(in) ::  & 
     34         tab2d 
     35 
     36      INTEGER, DIMENSION ( ndim1d), INTENT ( in) :: & 
     37         tab_ind 
     38 
     39      REAL(wp), DIMENSION(ndim1d), INTENT ( out) ::  &  
     40         tab1d 
     41 
     42      INTEGER ::  & 
     43         jn , jid, jjd 
     44         
     45      DO jn = 1, ndim1d 
     46         jid        = MOD( tab_ind(jn) - 1, ndim2d_x ) + 1 
     47         jjd        = ( tab_ind(jn) - 1 ) / ndim2d_x + 1 
     48         tab1d( jn) = tab2d( jid, jjd) 
    4349      END DO  
    44       ! 
     50 
    4551   END SUBROUTINE tab_2d_1d_2 
    4652 
    4753 
    48    SUBROUTINE tab_1d_2d_2( kdim1d, ptab2d, ptab_ind, ptab1d, kdim2d_i, kdim2d_j ) 
    49       !!------------------------------------------------------------------- 
    50       INTEGER , INTENT(in   )                                ::   kdim1d 
    51       INTEGER , INTENT(in   )                                ::   kdim2d_i, kdim2d_j 
    52       INTEGER , INTENT(in   ), DIMENSION(kdim1d)             ::   ptab_ind 
    53       REAL(wp), INTENT(in   ), DIMENSION(kdim1d)             ::   ptab1d 
    54       REAL(wp), INTENT(  out), DIMENSION(kdim2d_i, kdim2d_j) ::   ptab2d 
    55       !! 
    56       INTEGER ::   jn, jid, jjd   ! dummy loop indices 
    57       !!------------------------------------------------------------------- 
    58       ! 
    59       DO jn = 1, kdim1d 
    60          jid = MOD( ptab_ind(jn) - 1, kdim2d_i) + 1 
    61          jjd =    ( ptab_ind(jn) - 1 ) / kdim2d_i  + 1 
    62          ptab2d(jid,jjd) = ptab1d(jn) 
     54   SUBROUTINE tab_1d_2d_2 ( ndim1d, tab2d, tab_ind, tab1d, ndim2d_x, ndim2d_y ) 
     55 
     56      INTEGER, INTENT ( in) :: & 
     57         ndim1d, ndim2d_x, ndim2d_y 
     58 
     59      INTEGER, DIMENSION (ndim1d) , INTENT (in) :: & 
     60         tab_ind 
     61 
     62      REAL(wp), DIMENSION(ndim1d), INTENT (in) ::  & 
     63         tab1d   
     64 
     65      REAL(wp), DIMENSION (ndim2d_x, ndim2d_y), INTENT ( out) :: & 
     66         tab2d 
     67 
     68      INTEGER :: & 
     69         jn, jid, jjd 
     70 
     71      DO jn = 1, ndim1d 
     72         jid             = MOD( tab_ind(jn) - 1, ndim2d_x) + 1 
     73         jjd             =    ( tab_ind(jn) - 1 ) / ndim2d_x  + 1 
     74         tab2d(jid, jjd) = tab1d( jn) 
    6375      END DO 
    64       ! 
     76 
    6577   END SUBROUTINE tab_1d_2d_2 
    6678 
    67 #else 
    68    !!---------------------------------------------------------------------- 
    69    !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    70    !!---------------------------------------------------------------------- 
    7179#endif 
    72  
    73    !!====================================================================== 
    7480END MODULE limtab_2 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90

    r1855 r1857  
    22   !!====================================================================== 
    33   !!                  ***  MODULE limthd_2   *** 
    4    !!   LIM 2 ice model : thermodynamics 
     4   !!              LIM thermo ice model : ice thermodynamic 
    55   !!====================================================================== 
    6    !! History :  LIM  ! 2000-01 (UCL) Original code 
     6   !! History :  1.0  ! 2000-01 (LIM) 
    77   !!            2.0  ! 2002-07 (C. Ethe, G. Madec) F90 
    8    !!             -   ! 2003-08 (C. Ethe)  add lim_thd_init 
    9    !!             -   ! 2008-08  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface 
     8   !!            2.0  ! 2003-08 (C. Ethe)  add lim_thd_init 
     9   !!             -   ! 2008-2008  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface 
    1010   !!--------------------------------------------------------------------- 
    1111#if defined key_lim2 
     
    1313   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_thd_2      : thermodynamics of sea ice 
     15   !!   lim_thd_2      : thermodynamic of sea ice 
    1616   !!   lim_thd_init_2 : initialisation of sea-ice thermodynamic 
    1717   !!---------------------------------------------------------------------- 
     
    2323   USE lib_mpp 
    2424   USE iom             ! IOM library 
    25    USE ice_2           ! LIM 2 sea-ice variables 
     25   USE ice_2           ! LIM sea-ice variables 
    2626   USE sbc_oce         !  
    2727   USE sbc_ice         !  
    28    USE dom_ice_2       ! LIM 2 sea-ice domain 
     28   USE thd_ice_2       ! LIM thermodynamic sea-ice variables 
     29   USE dom_ice_2       ! LIM sea-ice domain 
    2930   USE limthd_zdf_2 
    3031   USE limthd_lac_2 
     
    4849#  include "domzgr_substitute.h90" 
    4950#  include "vectopt_loop_substitute.h90" 
    50    !!---------------------------------------------------------------------- 
    51    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     51   !!-------- ------------------------------------------------------------- 
     52   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009)  
    5253   !! $Id$ 
    5354   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    565566   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    566567   !!---------------------------------------------------------------------- 
     568CONTAINS 
     569   SUBROUTINE lim_thd_2         ! Dummy routine 
     570   END SUBROUTINE lim_thd_2 
    567571#endif 
    568572 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_lac_2.F90

    r1855 r1857  
    11MODULE limthd_lac_2 
     2#if defined key_lim2 
    23   !!====================================================================== 
    34   !!                       ***  MODULE limthd_lac_2   *** 
    4    !!   LIM 2 ice model :   thermodynamics -- lateral accretion of the ice  
    5    !!====================================================================== 
    6    !! History :  LIM  !  2001-04 (UCL)  original code 
    7    !!            2.0  !  2002-08 (C. Ethe, G. Madec)  F90, mpp 
     5   !!                lateral thermodynamic growth of the ice  
     6   !!====================================================================== 
     7 
    88   !!---------------------------------------------------------------------- 
    9 #if defined key_lim2 
    10    !!---------------------------------------------------------------------- 
    11    !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    12    !!---------------------------------------------------------------------- 
    13    !!   lim_lat_acr_2   : lateral accretion of ice 
    14    !!---------------------------------------------------------------------- 
     9   !!   lim_lat_acr_2    : lateral accretion of ice 
     10   !! * Modules used 
    1511   USE par_oce          ! ocean parameters 
    16    USE phycst           ! physical constants 
    17    USE ice_2            ! LIM 2 sea-ice variables 
    18    USE limistate_2      ! LIM 2 initial state 
     12   USE phycst 
     13   USE thd_ice_2 
     14   USE ice_2 
     15   USE limistate_2  
    1916      
    2017   IMPLICIT NONE 
    2118   PRIVATE 
    2219 
    23    PUBLIC   lim_thd_lac_2   ! called by lim_thd_2 
    24  
    25    REAL(wp) ::   epsi20 = 1.e-20   ! constant values 
    26    REAL(wp) ::   epsi13 = 1.e-13   ! 
    27    REAL(wp) ::   rzero  = 0.e0     ! 
    28    REAL(wp) ::   rone   = 1.e0     ! 
    29  
     20   !! * Routine accessibility 
     21   PUBLIC lim_thd_lac_2   ! called by lim_thd_2 
     22 
     23   !! * Module variables 
     24   REAL(wp)  ::           &  ! constant values 
     25      epsi20 = 1.e-20  ,  & 
     26      epsi13 = 1.e-13  ,  & 
     27      zzero  = 0.e0    ,  & 
     28      zone   = 1.e0 
    3029   !!---------------------------------------------------------------------- 
    31    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     30   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    3231   !! $Id$ 
    33    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     32   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    3433   !!---------------------------------------------------------------------- 
    3534CONTAINS 
     
    6160      !!               update h_snow_1d, h_ice_1d and tbif_1d(:,:)       
    6261      !!  
    63       !! References :   M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid 
    64       !!                Fichefet T. and M. Maqueda 1997, J. Geo. Res., 102(C6), 12609 -12646    
     62      !! ** References : 
     63      !!      M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid 
     64      !!      Fichefet T. and M. Maqueda 1997, J. Geo. Res., 102(C6),  
     65      !!                                                12609 -12646    
     66      !! History : 
     67      !!   1.0  !  01-04 (LIM)  original code 
     68      !!   2.0  !  02-08 (C. Ethe, G. Madec)  F90, mpp 
    6569      !!------------------------------------------------------------------- 
    66       INTEGER, INTENT(in) ::   kideb   ! start point on which the the computation is applied 
    67       INTEGER, INTENT(in) ::   kiut    ! end point on which the the computation is applied 
    68       !! 
    69       INTEGER ::   ji         !  dummy loop indices 
    70       INTEGER ::   iicefr     !  1 = existing ice ; 0 = no ice 
    71       INTEGER ::   iiceform   !  1 = ice formed   ; 0 = no ice formed 
    72       INTEGER ::   ihemis     !  hemisphere indicator  
    73       REAL(wp), DIMENSION(jpij) ::   zqbgow      !  heat budget of the open water (negative) 
    74       REAL(wp), DIMENSION(jpij) ::   zfrl_old    !  previous sea/ice fraction 
    75       REAL(wp), DIMENSION(jpij) ::   zhice_old   !  previous ice thickness 
    76       REAL(wp), DIMENSION(jpij) ::   zhice0      !  thickness of newly formed ice in leads 
    77       REAL(wp), DIMENSION(jpij) ::   zfrlmin     !  minimum fraction for leads 
    78       REAL(wp), DIMENSION(jpij) ::   zdhicbot    !  part of thickness of newly formed ice in leads which  
    79       !                                          !  has been already used in transport for example 
    80       REAL(wp) ::   zhemis     !  hemisphere (0 = North, 1 = South) 
    81       REAL(wp) ::   zhicenew   !  new ice thickness 
    82       REAL(wp) ::   zholds2    !  ratio of previous ice thickness and 2  
    83       REAL(wp) ::   zhnews2    !  ratio of new ice thickness and 2  
    84       REAL(wp) ::   zfrlnew    !  new sea/ice fraction 
    85       REAL(wp) ::   zfrld      !  ratio of sea/ice fraction and minimum fraction for leads 
    86       REAL(wp) ::   zfrrate    !  leads-closure rate 
    87       REAL(wp) ::   zdfrl      !  sea-ice fraction increment 
    88       REAL(wp) ::   zdh1 , zdh2, zdh3, zdh4, zdh5   ! tempory scalars 
    89       REAL(wp) ::   ztint, zta1, zta2, zta3, zta4   !    -       - 
    90       REAL(wp) ::   zah, zalpha, zbeta              !    -       - 
     70      !! * Arguments 
     71      INTEGER , INTENT(IN)::  & 
     72         kideb          ,   &  ! start point on which the the computation is applied 
     73         kiut                  ! end point on which the the computation is applied 
     74 
     75      !! * Local variables 
     76      INTEGER ::            & 
     77         ji             ,   &  !  dummy loop indices 
     78         iicefr         ,   &  !  1 = existing ice ; 0 = no ice 
     79         iiceform       ,   &  !  1 = ice formed   ; 0 = no ice formed 
     80         ihemis                !  dummy indice 
     81      REAL(wp), DIMENSION(jpij) :: & 
     82         zqbgow           ,  &  !  heat budget of the open water (negative) 
     83         zfrl_old         ,  &  !  previous sea/ice fraction 
     84         zhice_old        ,  &  !  previous ice thickness 
     85         zhice0           ,  &  !  thickness of newly formed ice in leads 
     86         zfrlmin          ,  &  !  minimum fraction for leads 
     87         zdhicbot               !  part of thickness of newly formed ice in leads which  
     88                                !  has been already used in transport for example 
     89      REAL(wp)  ::  & 
     90         zhemis           ,  &  !  hemisphere (0 = North, 1 = South) 
     91         zhicenew         ,  &  !  new ice thickness 
     92         zholds2          ,  &  !  ratio of previous ice thickness and 2  
     93         zhnews2          ,  &  !  ratio of new ice thickness and 2  
     94         zfrlnew          ,  &  !  new sea/ice fraction 
     95         zfrld            ,  &  !  ratio of sea/ice fraction and minimum fraction for leads 
     96         zfrrate          ,  &  !  leads-closure rate 
     97         zdfrl                  !  sea-ice fraction increment 
     98      REAL(wp)  ::  & 
     99         zdh1 , zdh2 , zdh3 , zdh4, zdh5   , &   ! tempory scalars 
     100         ztint , zta1 , zta2 , zta3 , zta4 , & 
     101         zah, zalpha , zbeta 
    91102      !!---------------------------------------------------------------------       
    92103                    
    93104      !-------------------------------------------------------------- 
    94105      !   Computation of the heat budget of the open water (negative) 
    95       !--------------------------------------------------------------  
     106      !-------------------------------------------------------------- 
     107       
    96108      DO ji = kideb , kiut       
    97109         zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) 
     
    119131      !     Morales Maqueda, 1995 - Fichefet and Morales Maqueda, 1997 
    120132      !--------------------------------------------------------------------- 
     133       
    121134!CDIR NOVERRCHK 
    122135      DO ji = kideb , kiut 
     
    132145         !--computation of the remaining part of ice thickness which has been already used 
    133146         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) ) &  
    134             &         -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic )  
     147                      -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic )  
    135148      END DO 
    136149  
     
    142155      !         (1-Anew) * hsnew = (1-Aold) * hsold             
    143156      !-------------------------------------------------------------------------------------------- 
     157       
    144158      DO ji = kideb , kiut 
    145159         iicefr       = 1 - MAX( 0, INT( SIGN( 1.5 * zone , zfrl_old(ji) - 1.0 + epsi13 ) ) ) 
     
    155169      !   Ajustement of ice internal temperatures 
    156170      !------------------------------------------------------- 
     171       
    157172      DO ji = kideb , kiut 
    158173         iicefr      = 1 - MAX( 0, INT( SIGN( 1.5 * zone , zfrl_old(ji) - 1.0 + epsi13 ) ) ) 
     
    198213      !           dV = Vnew - Vold 
    199214      !------------------------------------------------------------- 
     215       
    200216      DO ji = kideb , kiut 
    201217         dvlbq_1d  (ji) = ( 1. - frld_1d(ji) ) * h_ice_1d(ji) - ( 1. - zfrl_old(ji) ) * zhice_old(ji) 
    202218         rdmicif_1d(ji) = rdmicif_1d(ji) + rhoic * dvlbq_1d(ji) 
    203219      END DO 
    204       ! 
     220       
    205221   END SUBROUTINE lim_thd_lac_2 
    206     
    207222#else 
    208    !!---------------------------------------------------------------------- 
    209    !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    210    !!---------------------------------------------------------------------- 
     223   !!====================================================================== 
     224   !!                       ***  MODULE limthd_lac_2   *** 
     225   !!                           no sea ice model 
     226   !!====================================================================== 
     227CONTAINS 
     228   SUBROUTINE lim_thd_lac_2           ! Empty routine 
     229   END SUBROUTINE lim_thd_lac_2 
    211230#endif 
    212  
    213    !!====================================================================== 
    214231END MODULE limthd_lac_2 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90

    r1855 r1857  
    44   !!                thermodynamic growth and decay of the ice  
    55   !!====================================================================== 
    6    !! History :  1.0  !  2001-04 (LIM) Original code 
    7    !!            2.0  !  2002-08 (C. Ethe, G. Madec) F90 
     6   !! History :  1.0  !  01-04 (LIM) Original code 
     7   !!            2.0  !  02-08 (C. Ethe, G. Madec) F90 
    88   !!---------------------------------------------------------------------- 
    99#if defined key_lim2 
     
    1111   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
     13   !!---------------------------------------------------------------------- 
    1314   !!   lim_thd_zdf_2 : vertical accr./abl. and lateral ablation of sea ice 
    1415   !!---------------------------------------------------------------------- 
     16   !! * Modules used 
    1517   USE par_oce          ! ocean parameters 
    16    USE phycst           ! physical constants 
    17    USE ice_2            ! LIM 2 sea-ice variables 
    18    USE limistate_2      ! LIM 2 initial state 
    19    USE in_out_manager   ! I/O manager 
    20    !! 
     18   USE phycst           ! ??? 
     19   USE thd_ice_2 
     20   USE ice_2 
     21   USE limistate_2 
     22   USE in_out_manager 
    2123   USE cpl_oasis3, ONLY : lk_cpl 
    2224       
     
    2426   PRIVATE 
    2527 
    26    PUBLIC   lim_thd_zdf_2   ! called by lim_thd_2 
    27  
    28    REAL(wp) ::   epsi20 = 1.e-20   ! constant values 
    29    REAL(wp) ::   epsi13 = 1.e-13   ! 
    30    REAL(wp) ::   rzero  = 0.e0     ! 
    31    REAL(wp) ::   rone   = 1.e0     ! 
    32  
    33    !!---------------------------------------------------------------------- 
    34    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     28   PUBLIC   lim_thd_zdf_2        ! called by lim_thd_2 
     29 
     30   REAL(wp) ::   epsi20 = 1.e-20  ,  &  ! constant values 
     31      &          epsi13 = 1.e-13  ,  & 
     32      &          zzero  = 0.e0    ,  & 
     33      &          zone   = 1.e0 
     34   !!---------------------------------------------------------------------- 
     35   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    3536   !! $Id$ 
    3637   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    7273      !! 
    7374      INTEGER ::   ji       ! dummy loop indices 
    74       REAL(wp), DIMENSION(jpij,2) ::   zqcmlt   ! energy due to surface( /1 ) and bottom melting( /2 ) 
    75       REAL(wp), DIMENSION(jpij) ::   ztsmlt     ! snow/ice surface melting temperature 
    76       REAL(wp), DIMENSION(jpij) ::   ztbif      ! int. temp. at the mid-point of the 1st layer of the snow/ice sys.  
    77       REAL(wp), DIMENSION(jpij) ::   zksn       ! effective conductivity of snow 
    78       REAL(wp), DIMENSION(jpij) ::   zkic       ! effective conductivity of ice 
    79       REAL(wp), DIMENSION(jpij) ::   zksndh     ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys.  
    80       REAL(wp), DIMENSION(jpij) ::   zfcsu      ! conductive heat flux at the surface of the snow/ice system  
    81       REAL(wp), DIMENSION(jpij) ::   zfcsudt    ! = zfcsu * dt 
    82       REAL(wp), DIMENSION(jpij) ::   zi0        ! frac. of the net SW rad. which is not absorbed at the surface 
    83       REAL(wp), DIMENSION(jpij) ::   z1mi0      ! fraction of the net SW radiation absorbed at the surface 
    84       REAL(wp), DIMENSION(jpij) ::   zqmax      ! maximum energy stored in brine pockets 
    85       REAL(wp), DIMENSION(jpij) ::   zrcpdt     ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer) 
    86       REAL(wp), DIMENSION(jpij) ::   zts_old    ! previous surface temperature 
    87       REAL(wp), DIMENSION(jpij) ::   zidsn , z1midsn , zidsnic ! tempory variables 
    88       !! 
    89       REAL(wp), DIMENSION(jpij) ::   zfnet       ! net heat flux at the top surface( incl. conductive heat flux) 
    90       REAL(wp), DIMENSION(jpij) ::   zsprecip    ! snow accumulation 
    91       REAL(wp), DIMENSION(jpij) ::   zhsnw_old   ! previous snow thickness 
    92       REAL(wp), DIMENSION(jpij) ::   zdhictop    ! change in ice thickness due to top surf ablation/accretion 
    93       REAL(wp), DIMENSION(jpij) ::   zdhicbot    ! change in ice thickness due to bottom surf abl/acc 
    94       REAL(wp), DIMENSION(jpij) ::   zqsup       ! energy transmitted to ocean (coming from top surf abl/acc) 
    95       REAL(wp), DIMENSION(jpij) ::   zqocea      ! energy transmitted to ocean (coming from bottom sur abl/acc) 
    96       REAL(wp), DIMENSION(jpij) ::   zfrl_old    ! previous sea/ice fraction 
    97       REAL(wp), DIMENSION(jpij) ::   zfrld_1d    ! new sea/ice fraction 
    98       REAL(wp), DIMENSION(jpij) ::   zep         ! internal temperature of the 2nd layer of the snow/ice system 
    99       !! 
    100       REAL(wp), DIMENSION(3) ::   zplediag       ! principle diagonal, subdiag. and supdiag. of the  
    101       REAL(wp), DIMENSION(3) ::   zsubdiag       ! tri-diagonal matrix coming from the computation 
    102       REAL(wp), DIMENSION(3) ::   zsupdiag       ! of the temperatures inside the snow-ice system 
    103       REAL(wp), DIMENSION(3) ::   zsmbr          ! second member 
    104       REAL(wp) ::   zhsu        ! thickness of surface layer 
    105       REAL(wp) ::   zhe         ! effective thickness for compu. of equ. thermal conductivity 
    106       REAL(wp) ::   zheshth     ! = zhe / thth 
    107       REAL(wp) ::   zghe        ! correction factor of the thermal conductivity 
    108       REAL(wp) ::   zumsb       ! parameter for numerical method to solve heat-diffusion eq. 
    109       REAL(wp) ::   zkhsn       ! conductivity at the snow layer 
    110       REAL(wp) ::   zkhic       ! conductivity at the ice layers 
    111       REAL(wp) ::   zkint       ! equivalent conductivity at the snow-ice interface 
    112       REAL(wp) ::   zkhsnint    ! = zkint*dt / (hsn*rhosn*cpsn)   
    113       REAL(wp) ::   zkhicint    ! = 2*zkint*dt / (hic*rhoic*cpic) 
    114       REAL(wp) ::   zpiv1 , zpiv2        ! tempory scalars used to solve the tri-diagonal system 
    115       REAL(wp) ::   zb2, zd2, zb3, zd3 
    116       REAL(wp) ::   ztint          ! equivalent temperature at the snow-ice interface 
    117       REAL(wp) ::   zexp         ! exponential function of the ice thickness 
    118       REAL(wp) ::   zfsab        ! part of solar radiation stored in brine pockets 
    119       REAL(wp) ::   zfts         ! value of energy balance function when the temp. equal surf. temp. 
    120       REAL(wp) ::   zdfts        ! value of derivative of ztfs when the temp. equal surf. temp. 
    121       REAL(wp) ::   zdts         ! surface temperature increment 
    122       REAL(wp) ::   zqsnw_mlt    ! energy needed to melt snow 
    123       REAL(wp) ::   zdhsmlt      ! change in snow thickness due to melt 
    124       REAL(wp) ::   zhsn         ! snow thickness (previous+accumulation-melt) 
    125       REAL(wp) ::   zqsn_mlt_rem   ! remaining heat coming from snow melting 
    126       REAL(wp) ::   zqice_top_mlt   ! energy used to melt ice at top surface 
    127       REAL(wp) ::   zdhssub        ! change in snow thick. due to sublimation or evaporation 
    128       REAL(wp) ::   zdhisub        ! change in ice thick. due to sublimation or evaporation     
    129       REAL(wp) ::   zdhsn          ! snow ice thickness increment 
    130       REAL(wp) ::   zdtsn          ! snow internal temp. increment 
    131       REAL(wp) ::   zdtic          ! ice internal temp. increment 
    132       REAL(wp) ::   zqnes          ! conductive energy due to ice melting in the first ice layer 
    133       !! 
    134       REAL(wp) ::   ztbot           ! temperature at the bottom surface 
    135       REAL(wp) ::   zfcbot          ! conductive heat flux at bottom surface 
    136       REAL(wp) ::   zqice_bot       ! energy used for bottom melting/growing 
    137       REAL(wp) ::   zqice_bot_mlt   ! energy used for bottom melting 
    138       REAL(wp) ::   zqstbif_bot    ! part of energy stored in brine pockets used for bottom melting 
    139       REAL(wp) ::   zqstbif_old    ! tempory var. for zqstbif_bot 
    140       REAL(wp) ::   zdhicmlt        ! change in ice thickness due to bottom melting 
    141       REAL(wp) ::   zdhicm          ! change in ice thickness var.  
    142       REAL(wp) ::   zdhsnm          ! change in snow thickness var.  
    143       REAL(wp) ::   zhsnfi          ! snow thickness var.  
    144       REAL(wp) ::   zc1, zpc1, zc2, zpc2, zp1, zp2   ! temporary scalars 
    145       REAL(wp) ::   ztb2, ztb3 
    146       !! 
    147       REAL(wp) ::   zdrmh           ! change in snow/ice thick. after snow-ice formation 
    148       REAL(wp) ::   zhicnew         ! new ice thickness 
    149       REAL(wp) ::   zhsnnew         ! new snow thickness 
    150       REAL(wp) ::   zquot , ztneq   ! tempory temp. variables 
    151       REAL(wp) ::   zqice, zqicetot & ! total heat inside the snow/ice system 
    152       REAL(wp) ::   zdfrl           ! change in ice concentration 
    153       REAL(wp) ::   zdvsnvol        ! change in snow volume 
    154       REAL(wp) ::   zdrfrl1, zdrfrl2   ! tempory scalars 
    155       REAL(wp) ::   zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind 
    156       !!---------------------------------------------------------------------- 
    157  
    158       !----------------------------------------------------------------------- 
    159       !  1. Boundaries conditions for snow/ice system internal temperature 
    160       !       - If tbif_1d(ji,1) > rt0_snow, tbif_1d(ji,1) = rt0_snow  
    161       !       - If tbif_1d(ji,2/3) > rt0_ice, tbif_1d(ji,2/3) = rt0_ice  
    162       !     Computation of energies due to surface and bottom melting  
    163       !----------------------------------------------------------------------- 
    164       DO ji = kideb , kiut 
    165          zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 
    166          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) ) 
    167          !--computation of energy due to surface melting 
    168          zqcmlt(ji,1) = ( MAX ( zzero ,  & 
    169             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 
    170          !--computation of energy due to bottom melting 
    171          zqcmlt(ji,2) = ( MAX( zzero , & 
    172             &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
    173             &           + MAX( zzero , & 
    174             &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
    175             &           ) * ( 1.0 - zihic  ) 
    176          !--limitation of  snow/ice system internal temperature 
    177          tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) ) 
    178          tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) ) 
    179          tbif_1d(ji,3)   = MIN( rt0_ice , tbif_1d(ji,3) ) 
    180       END DO 
    181  
    182       !------------------------------------------- 
    183       !  2. Calculate some intermediate variables.   
    184       !------------------------------------------- 
    185       zhsu = hnzst         ! initialisation of the thickness of surface layer 
    186       ! 
    187       DO ji = kideb , kiut 
    188          zind   = MAX(  zzero , SIGN( zone , zhsu   - h_snow_1d(ji) )  ) 
    189          zihsn  = MAX(  zzero , SIGN( zone , hsndif - h_snow_1d(ji) )  ) 
    190          zihsn  = MAX(  zihsn , zind                                   ) 
    191          zihic  = MAX(  zzero , sign( zone , hicdif - h_ice_1d (ji) )  ) 
    192          ! 
    193          !     2.1. Computation of surface melting temperature 
    194          !---------------------------------------------------- 
    195          zind  = MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) ) 
    196          ztsmlt(ji) = ( 1.0 - zind ) * rt0_snow + zind * rt0_ice 
    197          ! 
    198          !     2.2. Effective conductivity of snow and ice 
    199          !----------------------------------------------- 
    200  
    201          !---computation of the correction factor on the thermal conductivity 
    202          !-- (Morales Maqueda, 1995 ; Fichefet and Morales Maqueda, 1997) 
    203          zhe      =  ( rcdsn / ( rcdsn + rcdic ) ) * h_ice_1d(ji)   & 
    204             &     + ( rcdic / ( rcdsn + rcdic ) ) * h_snow_1d(ji)  
    205          zihe     = MAX( zzero , SIGN( zone , 2.0 * zhe - thth ) ) 
    206          zheshth  = zhe / thth 
    207          zghe     = ( 1.0 - zihe ) * zheshth * ( 2.0 - zheshth )   & 
    208             &     +         zihe   * 0.5 * ( 1.5 + LOG( 2.0 * zheshth ) ) 
    209  
    210          !---effective conductivities  
    211          zksn(ji)  = zghe * rcdsn   
    212          zkic(ji)  = zghe * rcdic 
    213          ! 
    214          !     2.3. Computation of the conductive heat flux from the snow/ice  
    215          !          system interior toward the top surface 
    216          !------------------------------------------------------------------ 
    217  
    218          !---Thermal conductivity at the mid-point of the first snow/ice system layer 
    219          zksndh(ji) =   ( ( 1.0 - zihsn ) * 2.0 * zksn(ji) + zihsn * 4.0 * zkic(ji) )   & 
    220             &         / ( ( 1.0 - zihsn ) *  h_snow_1d(ji)                              & 
    221             &           +        zihsn   *  ( ( 1.0 + 3.0 * zihic ) * h_ice_1d(ji)      & 
    222             &           + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) ) 
    223  
    224          !---internal temperature at the mid-point of the first snow/ice system layer 
    225          ztbif(ji)  = ( 1.0 - zihsn ) * tbif_1d(ji,1)                       & 
    226             &       +         zihsn   * ( ( 1.0 - zihic ) * tbif_1d(ji,2)   & 
    227             &       +         zihic   * tfu_1d(ji)   ) 
    228          !---conductive heat flux  
    229          zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 
    230  
    231       END DO 
    232  
    233       !-------------------------------------------------------------------- 
    234       !  3. Calculate :  
    235       !     - fstbif_1d, part of solar radiation absorbing inside the ice 
    236       !       assuming an exponential absorption (Grenfell and Maykut, 1977) 
    237       !     - zqmax,  maximum energy stored in brine pockets 
    238       !     - qstbif_1d, total energy stored in brine pockets (updating) 
    239       !------------------------------------------------------------------- 
    240  
    241       DO ji = kideb , kiut 
    242          zihsn  = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) ) 
    243          zihic  = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) )      
    244          zind   = MAX( zzero , SIGN (zone , hicdif - h_ice_1d(ji) ) ) 
    245          !--Computation of the fraction of the net shortwave radiation which 
    246          !--penetrates inside the ice cover ( See Forcat) 
    247          zi0(ji)  = zihsn * ( fr1_i0_1d(ji) + zihic * fr2_i0_1d(ji) ) 
    248          zexp     = MIN( zone , EXP( -1.5 * ( h_ice_1d(ji) - zhsu ) ) ) 
    249          fstbif_1d(ji) = zi0(ji) * qsr_ice_1d(ji) * zexp 
    250          !--Computation of maximum energy stored in brine pockets zqmax and update 
    251          !--the total energy stored in brine pockets, if less than zqmax 
    252          zqmax(ji) = MAX( zzero , 0.5 * xlic * ( h_ice_1d(ji) - hicmin ) ) 
    253          zfsab   = zi0(ji) * qsr_ice_1d(ji) * ( 1.0 - zexp ) 
    254          zihq    = ( 1.0 - zind ) * MAX(zzero, SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) ) & 
    255             &    +         zind   * zone 
    256          qstbif_1d(ji) = ( qstbif_1d(ji) + ( 1.0 - zihq ) * zfsab * rdt_ice ) * swiqst 
    257          !--fraction of shortwave radiation absorbed at surface 
    258          ziexp = zihq * zexp + ( 1.0 - zihq ) * ( swiqst + ( 1.0 - swiqst ) * zexp ) 
    259          z1mi0(ji) = 1.0 - zi0(ji) * ziexp 
    260       END DO 
    261  
    262       !-------------------------------------------------------------------------------- 
    263       !  4. Computation of the surface temperature : determined by considering the  
    264       !     budget of a thin layer of thick. zhsu at the top surface (H. Grenier, 1995)  
    265       !     and based on a surface energy balance :  
    266       !     hsu * rcp * dT/dt = Fsr + Fnsr(T) + Fcs(T), 
    267       !     where - Fsr is the net absorbed solar radiation,  
    268       !           - Fnsr is the total non solar radiation (incoming and outgoing long-wave, 
    269       !             sensible and latent heat fluxes) 
    270       !           - Fcs the conductive heat flux at the top of surface 
    271       !------------------------------------------------------------------------------ 
    272  
    273       !     4.1. Computation of intermediate values 
    274       !--------------------------------------------- 
    275       DO ji = kideb, kiut 
    276          zrcpdt(ji) = ( rcpsn * MIN( h_snow_1d(ji) , zhsu )    & 
    277             &       + rcpic * MAX( zhsu - h_snow_1d(ji) , zzero ) ) / rdt_ice 
    278          zts_old(ji) =  sist_1d(ji) 
    279       END DO 
    280  
    281       !     4.2. Computation of surface temperature by expanding the eq. of energy balance 
    282       !          with Ts = Tp + DT. One obtain , F(Tp) + DT * DF(Tp) = 0 
    283       !          where  - F(Tp) = Fsr + Fnsr(Tp) + Fcs(Tp)  
    284       !                 - DF(Tp)= (dFnsr(Tp)/dT) + (dFcs(Tp)/dT) - hsu*rcp/dt 
    285       !--------------------------------------------------------------------------------- 
    286  
    287       DO ji = kideb, kiut 
    288          !---computation of the derivative of energy balance function  
    289          zdfts    =  zksndh(ji)   & ! contribution of the conductive heat flux 
    290             &      + zrcpdt(ji)   & ! contribution of hsu * rcp / dt 
    291             &      - dqns_ice_1d (ji)     ! contribution of the total non solar radiation  
    292          !---computation of the energy balance function  
    293          zfts    = - z1mi0 (ji) * qsr_ice_1d(ji)   & ! net absorbed solar radiation 
    294             &      - qns_ice_1d(ji)                & ! total non solar radiation 
    295             &      - zfcsu (ji)                      ! conductive heat flux from the surface 
    296          !---computation of surface temperature increment   
    297          zdts    = -zfts / zdfts 
    298          !---computation of the new surface temperature  
    299          sist_1d(ji) = sist_1d(ji) + zdts 
    300       END DO 
    301  
    302       !---------------------------------------------------------------------------- 
    303       !  5. Boundary condition at the top surface 
    304       !--    IF Tsb < Tmelt, Fnet = Fcs (the net heat flux equal the conductive heat flux) 
    305       !      Otherwise Tsb = Tmelt and Qnet(Tmelt) > 0  
    306       !      Fnet(Tmelt) is therefore the net surface flux needed for melting 
    307       !---------------------------------------------------------------------------- 
     75      REAL(wp), DIMENSION(jpij,2) ::   zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 ) 
     76      REAL(wp), DIMENSION(jpij) ::  & 
     77         ztsmlt      &    ! snow/ice surface melting temperature 
     78         ,ztbif      &    ! int. temp. at the mid-point of the 1st layer of the snow/ice sys.  
     79         ,zksn       &    ! effective conductivity of snow 
     80         ,zkic       &    ! effective conductivity of ice 
     81         ,zksndh     &    ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys.  
     82         , zfcsu     &    ! conductive heat flux at the surface of the snow/ice system  
     83         , zfcsudt   &    ! = zfcsu * dt 
     84         , zi0       &    ! frac. of the net SW rad. which is not absorbed at the surface 
     85         , z1mi0     &    ! fraction of the net SW radiation absorbed at the surface 
     86         , zqmax     &    ! maximum energy stored in brine pockets 
     87         , zrcpdt    &    ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer) 
     88         , zts_old   &    ! previous surface temperature 
     89         , zidsn , z1midsn , zidsnic ! tempory variables 
     90      REAL(wp), DIMENSION(jpij) ::   & 
     91          zfnet       &  ! net heat flux at the top surface( incl. conductive heat flux) 
     92          , zsprecip  &    ! snow accumulation 
     93          , zhsnw_old &    ! previous snow thickness 
     94          , zdhictop  &    ! change in ice thickness due to top surf ablation/accretion 
     95          , zdhicbot  &    ! change in ice thickness due to bottom surf abl/acc 
     96          , zqsup     &    ! energy transmitted to ocean (coming from top surf abl/acc) 
     97          , zqocea    &    ! energy transmitted to ocean (coming from bottom sur abl/acc) 
     98          , zfrl_old  &    ! previous sea/ice fraction 
     99          , zfrld_1d    &    ! new sea/ice fraction 
     100          , zep            ! internal temperature of the 2nd layer of the snow/ice system 
     101       REAL(wp), DIMENSION(3) :: &  
     102          zplediag  &    ! principle diagonal, subdiag. and supdiag. of the  
     103          , zsubdiag  &    ! tri-diagonal matrix coming from the computation 
     104          , zsupdiag  &    ! of the temperatures inside the snow-ice system 
     105          , zsmbr          ! second member 
     106       REAL(wp) :: &  
     107          zhsu     &     ! thickness of surface layer 
     108          , zhe      &     ! effective thickness for compu. of equ. thermal conductivity 
     109          , zheshth  &     ! = zhe / thth 
     110          , zghe     &     ! correction factor of the thermal conductivity 
     111          , zumsb    &     ! parameter for numerical method to solve heat-diffusion eq. 
     112          , zkhsn    &     ! conductivity at the snow layer 
     113          , zkhic    &     ! conductivity at the ice layers 
     114          , zkint    &     ! equivalent conductivity at the snow-ice interface 
     115          , zkhsnint &     ! = zkint*dt / (hsn*rhosn*cpsn)   
     116          , zkhicint &     ! = 2*zkint*dt / (hic*rhoic*cpic) 
     117          , zpiv1 , zpiv2  &       ! tempory scalars used to solve the tri-diagonal system 
     118          , zb2 , zd2 , zb3 , zd3 & 
     119          , ztint          ! equivalent temperature at the snow-ice interface 
     120       REAL(wp) :: &  
     121          zexp      &     ! exponential function of the ice thickness 
     122          , zfsab     &     ! part of solar radiation stored in brine pockets 
     123          , zfts      &     ! value of energy balance function when the temp. equal surf. temp. 
     124          , zdfts     &     ! value of derivative of ztfs when the temp. equal surf. temp. 
     125          , zdts      &     ! surface temperature increment 
     126          , zqsnw_mlt &     ! energy needed to melt snow 
     127          , zdhsmlt   &     ! change in snow thickness due to melt 
     128          , zhsn      &     ! snow thickness (previous+accumulation-melt) 
     129          , zqsn_mlt_rem &  ! remaining heat coming from snow melting 
     130          , zqice_top_mlt & ! energy used to melt ice at top surface 
     131          , zdhssub      &  ! change in snow thick. due to sublimation or evaporation 
     132          , zdhisub      &  ! change in ice thick. due to sublimation or evaporation     
     133          , zdhsn        &  ! snow ice thickness increment 
     134          , zdtsn        &  ! snow internal temp. increment 
     135          , zdtic        &  ! ice internal temp. increment 
     136          , zqnes          ! conductive energy due to ice melting in the first ice layer 
     137       REAL(wp) :: &  
     138          ztbot     &      ! temperature at the bottom surface 
     139          , zfcbot    &      ! conductive heat flux at bottom surface 
     140          , zqice_bot &      ! energy used for bottom melting/growing 
     141          , zqice_bot_mlt &  ! energy used for bottom melting 
     142          , zqstbif_bot  &  ! part of energy stored in brine pockets used for bottom melting 
     143          , zqstbif_old  &  ! tempory var. for zqstbif_bot 
     144          , zdhicmlt      &  ! change in ice thickness due to bottom melting 
     145          , zdhicm        &  ! change in ice thickness var.  
     146          , zdhsnm        &  ! change in snow thickness var.  
     147          , zhsnfi        &  ! snow thickness var.  
     148          , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables 
     149          , ztb2, ztb3 
     150       REAL(wp) :: &  
     151          zdrmh         &   ! change in snow/ice thick. after snow-ice formation 
     152          , zhicnew       &   ! new ice thickness 
     153          , zhsnnew       &   ! new snow thickness 
     154          , zquot , ztneq &   ! tempory temp. variables 
     155          , zqice, zqicetot & ! total heat inside the snow/ice system 
     156          , zdfrl         &   ! change in ice concentration 
     157          , zdvsnvol      &   ! change in snow volume 
     158          , zdrfrl1, zdrfrl2 &  ! tempory scalars 
     159          , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind 
     160       !!---------------------------------------------------------------------- 
     161 
     162       !----------------------------------------------------------------------- 
     163       !  1. Boundaries conditions for snow/ice system internal temperature 
     164       !       - If tbif_1d(ji,1) > rt0_snow, tbif_1d(ji,1) = rt0_snow  
     165       !       - If tbif_1d(ji,2/3) > rt0_ice, tbif_1d(ji,2/3) = rt0_ice  
     166       !     Computation of energies due to surface and bottom melting  
     167       !----------------------------------------------------------------------- 
     168        
     169       DO ji = kideb , kiut 
     170          zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 
     171          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) ) 
     172          !--computation of energy due to surface melting 
     173          zqcmlt(ji,1) = ( MAX ( zzero ,  & 
     174             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 
     175          !--computation of energy due to bottom melting 
     176          zqcmlt(ji,2) = ( MAX( zzero , & 
     177             &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
     178             &           + MAX( zzero , & 
     179             &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
     180             &           ) * ( 1.0 - zihic  ) 
     181          !--limitation of  snow/ice system internal temperature 
     182          tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) ) 
     183          tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) ) 
     184          tbif_1d(ji,3)   = MIN( rt0_ice , tbif_1d(ji,3) ) 
     185       END DO 
     186 
     187       !------------------------------------------- 
     188       !  2. Calculate some intermediate variables.   
     189       !------------------------------------------- 
     190        
     191       ! initialisation of the thickness of surface layer 
     192       zhsu = hnzst   
     193 
     194       DO ji = kideb , kiut 
     195          zind   = MAX( zzero , SIGN( zone , zhsu - h_snow_1d(ji) ) ) 
     196          zihsn  = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 
     197          zihsn  = MAX( zihsn , zind ) 
     198          zihic  = MAX( zzero , sign( zone , hicdif - h_ice_1d(ji) ) ) 
     199          !     2.1. Computation of surface melting temperature 
     200          !---------------------------------------------------- 
     201          zind  = MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) ) 
     202          ztsmlt(ji) = ( 1.0 - zind ) * rt0_snow + zind * rt0_ice 
     203          ! 
     204          !     2.2. Effective conductivity of snow and ice 
     205          !----------------------------------------------- 
     206 
     207          !---computation of the correction factor on the thermal conductivity 
     208          !-- (Morales Maqueda, 1995 ; Fichefet and Morales Maqueda, 1997) 
     209          zhe      =  ( rcdsn / ( rcdsn + rcdic ) ) * h_ice_1d(ji)   & 
     210             &     + ( rcdic / ( rcdsn + rcdic ) ) * h_snow_1d(ji)  
     211          zihe     = MAX( zzero , SIGN( zone , 2.0 * zhe - thth ) ) 
     212          zheshth  = zhe / thth 
     213          zghe     = ( 1.0 - zihe ) * zheshth * ( 2.0 - zheshth )   & 
     214             &     +         zihe   * 0.5 * ( 1.5 + LOG( 2.0 * zheshth ) ) 
     215 
     216          !---effective conductivities  
     217          zksn(ji)  = zghe * rcdsn   
     218          zkic(ji)  = zghe * rcdic 
     219 
     220          ! 
     221          !     2.3. Computation of the conductive heat flux from the snow/ice  
     222          !          system interior toward the top surface 
     223          !------------------------------------------------------------------ 
     224 
     225          !---Thermal conductivity at the mid-point of the first snow/ice system layer 
     226          zksndh(ji) =   ( ( 1.0 - zihsn ) * 2.0 * zksn(ji) + zihsn * 4.0 * zkic(ji) )   & 
     227             &         / ( ( 1.0 - zihsn ) *  h_snow_1d(ji)                              & 
     228             &           +        zihsn   *  ( ( 1.0 + 3.0 * zihic ) * h_ice_1d(ji)      & 
     229             &           + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) ) 
     230 
     231          !---internal temperature at the mid-point of the first snow/ice system layer 
     232          ztbif(ji)  = ( 1.0 - zihsn ) * tbif_1d(ji,1)                       & 
     233             &       +         zihsn   * ( ( 1.0 - zihic ) * tbif_1d(ji,2)   & 
     234             &       +         zihic   * tfu_1d(ji)   ) 
     235          !---conductive heat flux  
     236          zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 
     237 
     238       END DO 
     239 
     240       !-------------------------------------------------------------------- 
     241       !  3. Calculate :  
     242       !     - fstbif_1d, part of solar radiation absorbing inside the ice 
     243       !       assuming an exponential absorption (Grenfell and Maykut, 1977) 
     244       !     - zqmax,  maximum energy stored in brine pockets 
     245       !     - qstbif_1d, total energy stored in brine pockets (updating) 
     246       !------------------------------------------------------------------- 
     247 
     248       DO ji = kideb , kiut 
     249          zihsn  = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) ) 
     250          zihic  = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) )      
     251          zind   = MAX( zzero , SIGN (zone , hicdif - h_ice_1d(ji) ) ) 
     252          !--Computation of the fraction of the net shortwave radiation which 
     253          !--penetrates inside the ice cover ( See Forcat) 
     254          zi0(ji)  = zihsn * ( fr1_i0_1d(ji) + zihic * fr2_i0_1d(ji) ) 
     255          zexp     = MIN( zone , EXP( -1.5 * ( h_ice_1d(ji) - zhsu ) ) ) 
     256          fstbif_1d(ji) = zi0(ji) * qsr_ice_1d(ji) * zexp 
     257          !--Computation of maximum energy stored in brine pockets zqmax and update 
     258          !--the total energy stored in brine pockets, if less than zqmax 
     259          zqmax(ji) = MAX( zzero , 0.5 * xlic * ( h_ice_1d(ji) - hicmin ) ) 
     260          zfsab   = zi0(ji) * qsr_ice_1d(ji) * ( 1.0 - zexp ) 
     261          zihq    = ( 1.0 - zind ) * MAX(zzero, SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) ) & 
     262             &    +         zind   * zone 
     263          qstbif_1d(ji) = ( qstbif_1d(ji) + ( 1.0 - zihq ) * zfsab * rdt_ice ) * swiqst 
     264          !--fraction of shortwave radiation absorbed at surface 
     265          ziexp = zihq * zexp + ( 1.0 - zihq ) * ( swiqst + ( 1.0 - swiqst ) * zexp ) 
     266          z1mi0(ji) = 1.0 - zi0(ji) * ziexp 
     267       END DO 
     268 
     269       !-------------------------------------------------------------------------------- 
     270       !  4. Computation of the surface temperature : determined by considering the  
     271       !     budget of a thin layer of thick. zhsu at the top surface (H. Grenier, 1995)  
     272       !     and based on a surface energy balance :  
     273       !     hsu * rcp * dT/dt = Fsr + Fnsr(T) + Fcs(T), 
     274       !     where - Fsr is the net absorbed solar radiation,  
     275       !           - Fnsr is the total non solar radiation (incoming and outgoing long-wave, 
     276       !             sensible and latent heat fluxes) 
     277       !           - Fcs the conductive heat flux at the top of surface 
     278       !------------------------------------------------------------------------------ 
     279 
     280       !     4.1. Computation of intermediate values 
     281       !--------------------------------------------- 
     282       DO ji = kideb, kiut 
     283          zrcpdt(ji) = ( rcpsn * MIN( h_snow_1d(ji) , zhsu )    & 
     284             &       + rcpic * MAX( zhsu - h_snow_1d(ji) , zzero ) ) / rdt_ice 
     285          zts_old(ji) =  sist_1d(ji) 
     286       END DO 
     287 
     288       !     4.2. Computation of surface temperature by expanding the eq. of energy balance 
     289       !          with Ts = Tp + DT. One obtain , F(Tp) + DT * DF(Tp) = 0 
     290       !          where  - F(Tp) = Fsr + Fnsr(Tp) + Fcs(Tp)  
     291       !                 - DF(Tp)= (dFnsr(Tp)/dT) + (dFcs(Tp)/dT) - hsu*rcp/dt 
     292       !--------------------------------------------------------------------------------- 
     293 
     294       DO ji = kideb, kiut 
     295          !---computation of the derivative of energy balance function  
     296          zdfts    =  zksndh(ji)   & ! contribution of the conductive heat flux 
     297             &      + zrcpdt(ji)   & ! contribution of hsu * rcp / dt 
     298             &      - dqns_ice_1d (ji)     ! contribution of the total non solar radiation  
     299          !---computation of the energy balance function  
     300          zfts    = - z1mi0 (ji) * qsr_ice_1d(ji)   & ! net absorbed solar radiation 
     301             &      - qns_ice_1d(ji)                & ! total non solar radiation 
     302             &      - zfcsu (ji)                      ! conductive heat flux from the surface 
     303          !---computation of surface temperature increment   
     304          zdts    = -zfts / zdfts 
     305          !---computation of the new surface temperature  
     306          sist_1d(ji) = sist_1d(ji) + zdts 
     307       END DO 
     308 
     309       !---------------------------------------------------------------------------- 
     310       !  5. Boundary condition at the top surface 
     311       !--    IF Tsb < Tmelt, Fnet = Fcs (the net heat flux equal the conductive heat flux) 
     312       !      Otherwise Tsb = Tmelt and Qnet(Tmelt) > 0  
     313       !      Fnet(Tmelt) is therefore the net surface flux needed for melting 
     314       !---------------------------------------------------------------------------- 
    308315        
    309316        
    310       !     5.1.  Limitation of surface temperature and update total non solar fluxes, 
    311       !          latent heat flux and conductive flux at the top surface  
    312       !----------------------------------------------------------------------   
     317       !     5.1.  Limitation of surface temperature and update total non solar fluxes, 
     318       !          latent heat flux and conductive flux at the top surface  
     319       !----------------------------------------------------------------------   
    313320                      
    314       IF ( .NOT. lk_cpl ) THEN   ! duplicate the loop for performances issues 
    315          DO ji = kideb, kiut 
    316             sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 
    317             qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
    318             qla_ice_1d(ji) = qla_ice_1d(ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
    319             zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 
    320          END DO 
    321       ELSE 
    322          DO ji = kideb, kiut 
    323             sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 
    324             zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 
    325          END DO 
    326       ENDIF 
    327  
    328       !     5.2. Calculate available heat for surface ablation.  
    329       !--------------------------------------------------------------------- 
    330  
    331       DO ji = kideb, kiut 
    332          zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)           
    333          zfnet(ji) = MAX( zzero , zfnet(ji) ) 
    334          zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) ) 
    335       END DO 
    336  
    337       !------------------------------------------------------------------------- 
    338       !  6. Calculate changes in internal temperature due to vertical diffusion    
    339       !     processes. The evolution of this temperature is governed by the one- 
    340       !     dimensionnal heat-diffusion equation.  
    341       !     Given the temperature tbif(1/2/3), at time m we solve a set 
    342       !     of finite difference equations to obtain new tempe. Each tempe is coupled 
    343       !     to the temp. immediatly above and below by heat conduction terms. Thus  
    344       !     we have a set of equations of the form A * T = B, where A is a tridiagonal 
    345       !     matrix, T a vector whose components are the unknown new temp. 
    346       !------------------------------------------------------------------------- 
     321       IF ( .NOT. lk_cpl ) THEN   ! duplicate the loop for performances issues 
     322          DO ji = kideb, kiut 
     323             sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 
     324             qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
     325             qla_ice_1d(ji) = qla_ice_1d(ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 
     326             zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 
     327          END DO 
     328       ELSE 
     329          DO ji = kideb, kiut 
     330             sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 
     331             zfcsu(ji)  = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 
     332          END DO 
     333       ENDIF 
     334 
     335       !     5.2. Calculate available heat for surface ablation.  
     336       !--------------------------------------------------------------------- 
     337 
     338       DO ji = kideb, kiut 
     339          zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)           
     340          zfnet(ji) = MAX( zzero , zfnet(ji) ) 
     341          zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) ) 
     342       END DO 
     343 
     344       !------------------------------------------------------------------------- 
     345       !  6. Calculate changes in internal temperature due to vertical diffusion    
     346       !     processes. The evolution of this temperature is governed by the one- 
     347       !     dimensionnal heat-diffusion equation.  
     348       !     Given the temperature tbif(1/2/3), at time m we solve a set 
     349       !     of finite difference equations to obtain new tempe. Each tempe is coupled 
     350       !     to the temp. immediatly above and below by heat conduction terms. Thus  
     351       !     we have a set of equations of the form A * T = B, where A is a tridiagonal 
     352       !     matrix, T a vector whose components are the unknown new temp. 
     353       !------------------------------------------------------------------------- 
    347354        
    348       !--parameter for the numerical methode use to solve the heat-diffusion equation 
    349       !- implicit, explicit or Crank-Nicholson 
    350       zumsb = 1.0 - sbeta   
    351       DO ji = kideb, kiut 
    352          zidsn(ji)   = MAX ( zzero, SIGN( zone, hsndif - h_snow_1d(ji) ) )  
    353          z1midsn(ji) = 1.0 - zidsn(ji) 
    354          zihic       = MAX ( zzero, SIGN( zone, hicdif - h_ice_1d(ji) ) )  
    355          zidsnic(ji) = zidsn(ji) *  zihic  
    356          zfcsudt(ji) = zfcsu(ji) * rdt_ice  
    357       END DO  
     355       !--parameter for the numerical methode use to solve the heat-diffusion equation 
     356       !- implicit, explicit or Crank-Nicholson 
     357       zumsb = 1.0 - sbeta   
     358       DO ji = kideb, kiut 
     359          zidsn(ji)   = MAX ( zzero, SIGN( zone, hsndif - h_snow_1d(ji) ) )  
     360          z1midsn(ji) = 1.0 - zidsn(ji) 
     361          zihic       = MAX ( zzero, SIGN( zone, hicdif - h_ice_1d(ji) ) )  
     362          zidsnic(ji) = zidsn(ji) *  zihic  
     363          zfcsudt(ji) = zfcsu(ji) * rdt_ice  
     364       END DO  
    358365    
    359       DO ji = kideb, kiut 
    360  
    361          !     6.1 Calculate intermediate variables. 
    362          !---------------------------------------- 
    363  
    364          !--conductivity at the snow surface 
    365          zkhsn = 2.0 * zksn(ji) * rdt_ice / rcpsn 
    366          !--conductivity at the ice surface 
    367          zkhic = 4.0 * zkic(ji) * rdt_ice / MAX( h_ice_1d(ji) * h_ice_1d(ji) * rcpic , epsi20 ) 
    368          !--conductivity at the snow/ice interface  
    369          zkint = 4.0 * zksn(ji) * zkic(ji)  & 
    370             &        / ( zksn(ji) * h_ice_1d(ji) + 2.0 * zkic(ji) * h_snow_1d(ji) * z1midsn(ji))  
    371          zkhsnint = zkint * rdt_ice / rcpsn 
    372          zkhicint = zkint * 2.0 * rdt_ice / MAX( h_ice_1d(ji) * rcpic , epsi20 ) 
    373            
    374          !     6.2. Fulfill the linear system matrix. 
    375          !----------------------------------------- 
     366       DO ji = kideb, kiut 
     367 
     368          !     6.1 Calculate intermediate variables. 
     369          !---------------------------------------- 
     370 
     371          !--conductivity at the snow surface 
     372          zkhsn = 2.0 * zksn(ji) * rdt_ice / rcpsn 
     373          !--conductivity at the ice surface 
     374          zkhic = 4.0 * zkic(ji) * rdt_ice / MAX( h_ice_1d(ji) * h_ice_1d(ji) * rcpic , epsi20 ) 
     375          !--conductivity at the snow/ice interface  
     376          zkint = 4.0 * zksn(ji) * zkic(ji)  & 
     377             &        / ( zksn(ji) * h_ice_1d(ji) + 2.0 * zkic(ji) * h_snow_1d(ji) * z1midsn(ji))  
     378          zkhsnint = zkint * rdt_ice / rcpsn 
     379          zkhicint = zkint * 2.0 * rdt_ice / MAX( h_ice_1d(ji) * rcpic , epsi20 ) 
     380           
     381          !     6.2. Fulfill the linear system matrix. 
     382          !----------------------------------------- 
    376383!$$$          zplediag(1) = 1 + sbeta * z1midsn(ji) * ( zkhsn + zkhsnint )        
    377          zplediag(1) =   zidsn(ji) + z1midsn(ji) * h_snow_1d(ji)   & 
    378             &          + sbeta * z1midsn(ji) * zkhsnint  
    379          zplediag(2) = 1 + sbeta * ( z1midsn(ji) * zkhicint + zkhic )  
    380          zplediag(3) = 1 + 3.0 * sbeta * zkhic    
    381  
    382          zsubdiag(1) =  0.e0               
    383          zsubdiag(2) = -1.e0 * z1midsn(ji) * sbeta * zkhicint 
    384          zsubdiag(3) = -1.e0 * sbeta * zkhic  
    385  
    386          zsupdiag(1) = -1.e0 * z1midsn(ji) * sbeta * zkhsnint  
    387          zsupdiag(2) = zsubdiag(3) 
    388          zsupdiag(3) =  0.e0 
    389            
    390          !     6.3. Fulfill the idependent term vector. 
    391          !------------------------------------------- 
     384          zplediag(1) =   zidsn(ji) + z1midsn(ji) * h_snow_1d(ji)   & 
     385             &          + sbeta * z1midsn(ji) * zkhsnint  
     386          zplediag(2) = 1 + sbeta * ( z1midsn(ji) * zkhicint + zkhic )  
     387          zplediag(3) = 1 + 3.0 * sbeta * zkhic    
     388 
     389          zsubdiag(1) =  0.e0               
     390          zsubdiag(2) = -1.e0 * z1midsn(ji) * sbeta * zkhicint 
     391          zsubdiag(3) = -1.e0 * sbeta * zkhic  
     392 
     393          zsupdiag(1) = -1.e0 * z1midsn(ji) * sbeta * zkhsnint  
     394          zsupdiag(2) = zsubdiag(3) 
     395          zsupdiag(3) =  0.e0 
     396           
     397          !     6.3. Fulfill the idependent term vector. 
     398          !------------------------------------------- 
    392399           
    393400!$$$          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *   & 
     
    395402!$$$             &         - zumsb * ( zkhsn * tbif_1d(ji,1) 
    396403!$$$             &                   + zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) ) ) 
    397          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *    & 
    398             &       ( h_snow_1d(ji) * tbif_1d(ji,1) - ( zfcsudt(ji) / rcpsn )  & 
    399             &       - zumsb * zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) ) 
    400  
    401          zsmbr(2) =  tbif_1d(ji,2)  & 
    402             &      - zidsn(ji) * ( 1.0 - zidsnic(ji) ) & 
    403             &        * ( zfcsudt(ji) / MAX( h_ice_1d(ji) * rcpic , epsi20 ) ) & 
    404             &      + zumsb * ( zkhicint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) & 
    405             &                   - zkhic * ( tbif_1d(ji,2) - tbif_1d(ji,3) )  ) 
    406  
    407          zsmbr(3) =  tbif_1d(ji,3)  & 
    408             &      + zkhic * ( 2.0 * tfu_1d(ji) & 
    409             &                + zumsb * ( tbif_1d(ji,2) - 3.0 * tbif_1d(ji,3) ) ) 
    410            
    411          !     6.4. Solve the system (Gauss elimination method). 
    412          !---------------------------------------------------- 
    413            
    414          zpiv1 = zsubdiag(2) / zplediag(1)  
    415          zb2   = zplediag(2) - zpiv1 * zsupdiag(1) 
    416          zd2   = zsmbr(2) - zpiv1 * zsmbr(1) 
    417  
    418          zpiv2 = zsubdiag(3) / zb2 
    419          zb3   = zplediag(3) - zpiv2 * zsupdiag(2) 
    420          zd3   = zsmbr(3) - zpiv2 * zd2 
    421  
    422          tbif_1d(ji,3) = zd3 / zb3 
    423          tbif_1d(ji,2) = ( zd2 - zsupdiag(2) * tbif_1d(ji,3) ) / zb2 
    424          tbif_1d(ji,1) = ( zsmbr(1) - zsupdiag(1) * tbif_1d(ji,2) ) / zplediag(1)             
    425  
    426          !--- taking into account the particular case of  zidsnic(ji) = 1 
    427          ztint =  (  zkic(ji) * h_snow_1d(ji) * tfu_1d (ji)    & 
    428             &      + zksn(ji) * h_ice_1d(ji) * sist_1d(ji) )   & 
    429             &   / ( zkic(ji) * h_snow_1d(ji) + zksn(ji) * h_ice_1d(ji) )  
    430  
    431          tbif_1d(ji,1) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,1)   & 
    432             &                + zidsnic(ji)   * ( ztint + sist_1d(ji) ) / 2.0 
    433          tbif_1d(ji,2) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,2)   & 
    434             &                + zidsnic(ji)   * ( 3.0 * ztint + tfu_1d(ji) ) / 4.0 
    435          tbif_1d(ji,3) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3)   & 
    436             &                + zidsnic(ji)   * ( ztint + 3.0 * tfu_1d(ji) ) / 4.0      
    437       END DO 
     404          zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) *    & 
     405             &       ( h_snow_1d(ji) * tbif_1d(ji,1) - ( zfcsudt(ji) / rcpsn )  & 
     406             &       - zumsb * zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) ) 
     407 
     408          zsmbr(2) =  tbif_1d(ji,2)  & 
     409             &      - zidsn(ji) * ( 1.0 - zidsnic(ji) ) & 
     410             &        * ( zfcsudt(ji) / MAX( h_ice_1d(ji) * rcpic , epsi20 ) ) & 
     411             &      + zumsb * ( zkhicint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) & 
     412             &                   - zkhic * ( tbif_1d(ji,2) - tbif_1d(ji,3) )  ) 
     413 
     414          zsmbr(3) =  tbif_1d(ji,3)  & 
     415             &      + zkhic * ( 2.0 * tfu_1d(ji) & 
     416             &                + zumsb * ( tbif_1d(ji,2) - 3.0 * tbif_1d(ji,3) ) ) 
     417           
     418          !     6.4. Solve the system (Gauss elimination method). 
     419          !---------------------------------------------------- 
     420           
     421          zpiv1 = zsubdiag(2) / zplediag(1)  
     422          zb2   = zplediag(2) - zpiv1 * zsupdiag(1) 
     423          zd2   = zsmbr(2) - zpiv1 * zsmbr(1) 
     424 
     425          zpiv2 = zsubdiag(3) / zb2 
     426          zb3   = zplediag(3) - zpiv2 * zsupdiag(2) 
     427          zd3   = zsmbr(3) - zpiv2 * zd2 
     428 
     429          tbif_1d(ji,3) = zd3 / zb3 
     430          tbif_1d(ji,2) = ( zd2 - zsupdiag(2) * tbif_1d(ji,3) ) / zb2 
     431          tbif_1d(ji,1) = ( zsmbr(1) - zsupdiag(1) * tbif_1d(ji,2) ) / zplediag(1)             
     432 
     433          !--- taking into account the particular case of  zidsnic(ji) = 1 
     434          ztint =  (  zkic(ji) * h_snow_1d(ji) * tfu_1d (ji)    & 
     435             &      + zksn(ji) * h_ice_1d(ji) * sist_1d(ji) )   & 
     436             &   / ( zkic(ji) * h_snow_1d(ji) + zksn(ji) * h_ice_1d(ji) )  
     437 
     438          tbif_1d(ji,1) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,1)   & 
     439             &                + zidsnic(ji)   * ( ztint + sist_1d(ji) ) / 2.0 
     440          tbif_1d(ji,2) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,2)   & 
     441             &                + zidsnic(ji)   * ( 3.0 * ztint + tfu_1d(ji) ) / 4.0 
     442          tbif_1d(ji,3) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3)   & 
     443             &                + zidsnic(ji)   * ( ztint + 3.0 * tfu_1d(ji) ) / 4.0      
     444       END DO 
    438445  
    439       !---------------------------------------------------------------------- 
    440       !  9. Take into account surface ablation and bottom accretion-ablation.| 
    441       !---------------------------------------------------------------------- 
     446       !---------------------------------------------------------------------- 
     447       !  9. Take into account surface ablation and bottom accretion-ablation.| 
     448       !---------------------------------------------------------------------- 
    442449        
    443       !---Snow accumulation in one thermodynamic time step 
    444       zsprecip(kideb:kiut) = sprecip_1d(kideb:kiut) * rdt_ice / rhosn 
    445  
    446  
    447       DO ji = kideb, kiut 
    448            
    449          !      9.1. Surface ablation and update of snow thickness and qstbif_1d 
    450          !-------------------------------------------------------------------- 
     450       !---Snow accumulation in one thermodynamic time step 
     451       zsprecip(kideb:kiut) = sprecip_1d(kideb:kiut) * rdt_ice / rhosn 
     452 
     453 
     454       DO ji = kideb, kiut 
     455           
     456          !      9.1. Surface ablation and update of snow thickness and qstbif_1d 
     457          !-------------------------------------------------------------------- 
    451458           
    452459          !-------------------------------------------------------------------------- 
     
    752759          qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji) 
    753760          frld_1d(ji)    = zfrld_1d(ji) 
    754          ! 
    755       END DO 
    756       !  
     761          ! 
     762       END DO 
     763       !  
     764    END SUBROUTINE lim_thd_zdf_2 
     765 
     766#else 
     767   !!---------------------------------------------------------------------- 
     768   !!   Default Option                                     NO sea-ice model 
     769   !!---------------------------------------------------------------------- 
     770CONTAINS 
     771   SUBROUTINE lim_thd_zdf_2          ! Empty routine 
    757772   END SUBROUTINE lim_thd_zdf_2 
    758  
    759 #else 
    760    !!---------------------------------------------------------------------- 
    761    !!   Default Option                                     NO sea-ice model 
    762    !!---------------------------------------------------------------------- 
    763773#endif 
    764774 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtrp_2.F90

    r1855 r1857  
    44   !! LIM 2.0 transport ice model : sea-ice advection/diffusion 
    55   !!====================================================================== 
    6    !! History :  LIM  ! 2000-01 (LIM)  Original code 
    7    !!            1.0  ! 2001-05 (G. Madec, R. Hordoir) opa norm 
    8    !!            2.0  ! 2004-01 (G. Madec, C. Ethe)  F90, mpp 
    9    !!---------------------------------------------------------------------- 
    106#if defined key_lim2 
    117   !!---------------------------------------------------------------------- 
     
    1511   !!   lim_trp_init_2 : initialization and namelist read 
    1612   !!---------------------------------------------------------------------- 
     13   !! * Modules used 
    1714   USE phycst 
    1815   USE dom_oce 
     
    2926   PRIVATE 
    3027 
    31    PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2 
    32  
    33    REAL(wp), PUBLIC  ::   bound  = 0.e0   !: boundary condit. (0.0 no-slip, 1.0 free-slip) 
    34  
    35    REAL(wp) ::   epsi06 = 1.e-06   ! constant values 
    36    REAL(wp) ::   epsi03 = 1.e-03   ! 
    37    REAL(wp) ::   epsi16 = 1.e-16   ! 
    38    REAL(wp) ::   rzero  = 0.e0     ! 
    39    REAL(wp) ::   rone   = 1.e0     ! 
     28   !! * Routine accessibility 
     29   PUBLIC lim_trp_2     ! called by sbc_ice_lim_2 
     30 
     31   !! * Shared module variables 
     32   REAL(wp), PUBLIC  ::   &  !: 
     33      bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip) 
     34 
     35   !! * Module variables 
     36   REAL(wp)  ::           &  ! constant values 
     37      epsi06 = 1.e-06  ,  & 
     38      epsi03 = 1.e-03  ,  & 
     39      epsi16 = 1.e-16  ,  & 
     40      rzero  = 0.e0    ,  & 
     41      rone   = 1.e0 
    4042 
    4143   !! * Substitution 
    4244#  include "vectopt_loop_substitute.h90" 
    4345   !!---------------------------------------------------------------------- 
    44    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     46   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    4547   !! $Id$ 
    46    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     48   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    4749   !!---------------------------------------------------------------------- 
    4850 
     
    6062      !! 
    6163      !! ** action : 
     64      !! 
     65      !! History : 
     66      !!   1.0  !  00-01 (LIM)  Original code 
     67      !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
     68      !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp 
    6269      !!--------------------------------------------------------------------- 
    6370      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    64       !! 
    65       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    66       INTEGER  ::   initad       ! number of sub-timestep for the advection 
    67       REAL(wp) ::   zindb , zacrith, zindsn , zignm , zvbord, zrtt, ztic1, zusnit 
    68       REAL(wp) ::   zindic, zusvosn, zusvoic, zindhe, zcfl  , ztsn, ztic2 
    69       REAL(wp), DIMENSION(jpi,jpj) ::   zui_u , zs0sn, zs0c0, zs0a , zsm      ! 2D workspace 
    70       REAL(wp), DIMENSION(jpi,jpj) ::   zvi_v , zs0st, zs0c1, zs0c2, zs0ice   !  -      - 
     71 
     72      INTEGER  ::   ji, jj, jk,   &  ! dummy loop indices 
     73         &          initad           ! number of sub-timestep for the advection 
     74 
     75      REAL(wp) ::  &                               
     76         zindb  ,  & 
     77         zacrith, & 
     78         zindsn , & 
     79         zindic , & 
     80         zusvosn, & 
     81         zusvoic, & 
     82         zignm  , & 
     83         zindhe , & 
     84         zvbord , & 
     85         zcfl   , & 
     86         zusnit , & 
     87         zrtt, ztsn, ztic1, ztic2 
     88 
     89      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace 
     90         zui_u , zvi_v , zsm   ,         & 
     91         zs0ice, zs0sn , zs0a  ,         & 
     92         zs0c0 , zs0c1 , zs0c2 ,         & 
     93         zs0st 
    7194      !--------------------------------------------------------------------- 
    7295 
     
    7598      zsm(:,:) = area(:,:) 
    7699       
    77       !                             !-------------------------------------! 
    78       IF( ln_limdyn ) THEN          !   Advection of sea ice properties   ! 
    79          !                          !-------------------------------------! 
     100      IF( ln_limdyn ) THEN 
     101         !-------------------------------------! 
     102         !   Advection of sea ice properties   ! 
     103         !-------------------------------------! 
    80104 
    81105         ! ice velocities at ocean U- and V-points (zui_u,zvi_v) 
     
    89113            END DO 
    90114         END DO 
    91          CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )         ! Lateral boundary conditions 
     115         ! Lateral boundary conditions on zui_u, zvi_v 
     116         CALL lbc_lnk( zui_u, 'U', -1. ) 
     117         CALL lbc_lnk( zvi_v, 'V', -1. ) 
    92118 
    93119         ! CFL test for stability 
     
    97123         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) ) 
    98124 
    99          IF( lk_mpp ) CALL mpp_max( zcfl ) 
    100  
    101          IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : CFL violation at the ',nday,'th day, cfl = ',zcfl 
     125         IF (lk_mpp ) CALL mpp_max(zcfl) 
     126 
     127         IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl 
    102128 
    103129         ! content of properties 
     
    118144         zusnit = 1.0 / REAL( initad )  
    119145          
    120          IF( MOD( nday , 2 ) == 0) THEN 
     146         IF ( MOD( nday , 2 ) == 0) THEN 
    121147            DO jk = 1,initad 
    122148               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice ) 
     
    202228         !   Up-dating and limitation of sea ice properties after transport   ! 
    203229         ! -------------------------------------------------------------------! 
     230 
     231         ! Up-dating and limitation of sea ice properties after transport. 
    204232         DO jj = 1, jpj 
     233!!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq! 
    205234            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH 
    206235            DO ji = 1, jpi 
    207                !                                                        ! Recover mean values over the grid squares. 
     236 
     237               ! Recover mean values over the grid squares. 
    208238               zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) ) 
    209239               zs0ice(ji,jj) = MAX( rzero, zs0ice(ji,jj)/area(ji,jj) ) 
     
    213243               zs0c2 (ji,jj) = MAX( rzero, zs0c2 (ji,jj)/area(ji,jj) ) 
    214244               zs0st (ji,jj) = MAX( rzero, zs0st (ji,jj)/area(ji,jj) ) 
    215                !                                                        ! Recover in situ values. 
     245 
     246               ! Recover in situ values. 
    216247               zindb         = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) ) 
    217248               zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) ) 
     
    231262               zrtt          = 173.15 * rone  
    232263               ztsn          =          zignm   * tbif(ji,jj,1)  & 
    233                   &           + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) )  
     264                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) )  
    234265               ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) ) 
    235266               ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) ) 
    236                ! 
     267  
    237268               tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)                
    238269               tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj) 
     
    241272            END DO 
    242273         END DO 
    243          ! 
     274          
    244275      ENDIF 
    245       ! 
     276       
    246277   END SUBROUTINE lim_trp_2 
    247278 
     
    257288      !! 
    258289      !! ** input   :   Namelist namicetrp 
     290      !! 
     291      !! history : 
     292      !!   2.0  !  03-08 (C. Ethe)  Original code 
    259293      !!------------------------------------------------------------------- 
    260294      NAMELIST/namicetrp/ bound 
    261295      !!------------------------------------------------------------------- 
    262       ! 
    263       REWIND ( numnam_ice )                  ! Read Namelist namicetrp 
     296 
     297      ! Read Namelist namicetrp 
     298      REWIND ( numnam_ice ) 
    264299      READ   ( numnam_ice  , namicetrp ) 
    265       IF(lwp) THEN                           ! control print 
     300      IF(lwp) THEN 
    266301         WRITE(numout,*) 
    267302         WRITE(numout,*) 'lim_trp_init_2 : Ice parameters for advection ' 
     
    269304         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound 
    270305      ENDIF 
    271       !  
     306             
    272307   END SUBROUTINE lim_trp_init_2 
    273308 
     
    276311   !!   Default option         Empty Module                No sea-ice model 
    277312   !!---------------------------------------------------------------------- 
     313CONTAINS 
     314   SUBROUTINE lim_trp_2        ! Empty routine 
     315   END SUBROUTINE lim_trp_2 
    278316#endif 
    279317 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limwri_2.F90

    r1855 r1857  
    1212   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!   lim_wri_2       : write of the diagnostics variables in ouput file  
    15    !!   lim_wri_init_2  : initialization and namelist read 
    16    !!   lim_wri_state_2 : write for initial state (output.init.nc if ninist=1) or/and in the abort file 
     14   !!---------------------------------------------------------------------- 
     15   !!   lim_wri_2      : write of the diagnostics variables in ouput file  
     16   !!   lim_wri_init_2 : initialization and namelist read 
     17   !!   lim_wri_state_2 : write for initial state or/and abandon: 
     18   !!                     > output.init.nc (if ninist = 1 in namelist) 
     19   !!                     > output.abort.nc 
    1720   !!---------------------------------------------------------------------- 
    1821   USE phycst 
     
    3942   INTEGER, PARAMETER                       ::   jpnoumax = 40   ! maximum number of variable for ice output 
    4043   INTEGER                                  ::   noumef          ! number of fields 
    41    REAL(wp)           , DIMENSION(jpnoumax) ::   cmulti         ! multiplicative constant 
    42    REAL(wp)           , DIMENSION(jpnoumax) ::   cadd            ! additive constant 
     44   REAL(wp)           , DIMENSION(jpnoumax) ::   cmulti ,     &  ! multiplicative constant 
     45      &                                          cadd            ! additive constant 
    4346   CHARACTER(len = 35), DIMENSION(jpnoumax) ::   titn            ! title of the field 
    4447   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   nam             ! name of the field 
     
    4952   INTEGER , DIMENSION( jpij ) ::   ndex51              ! ???? 
    5053 
    51    REAL(wp) ::   epsi16 = 1.e-16   ! constant values 
    52    REAL(wp) ::   rzero  = 0.e0     ! 
    53    REAL(wp) ::   rone   = 1.e0     ! 
     54   REAL(wp)  ::            &  ! constant values 
     55      epsi16 = 1.e-16   ,  & 
     56      zzero  = 0.e0     ,  & 
     57      zone   = 1.e0 
    5458 
    5559   !! * Substitutions 
    5660#   include "vectopt_loop_substitute.h90" 
    5761   !!---------------------------------------------------------------------- 
    58    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     62   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006) 
    5963   !! $Id$ 
    6064   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    7781      !! 
    7882      !! ** Method  :   computes the average of some variables and write 
    79       !!              it in the NetCDF ouput files 
    80       !! CAUTION: the sea-ice time-step must be an integer fraction of a day 
     83      !!      it in the NetCDF ouput files 
     84      !!      CAUTION: the sea-ice time-step must be an integer fraction 
     85      !!      of a day 
    8186      !!------------------------------------------------------------------- 
    8287      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    8388      !! 
    84       INTEGER  ::   ji, jj, jf   ! dummy loop indices 
     89      INTEGER  ::   ji, jj, jf                      ! dummy loop indices 
    8590      CHARACTER(len = 40)  ::   clhstnam, clop 
    86       REAL(wp) ::   zsto, zjulian, zout          ! temporary scalars 
    87       REAL(wp) ::   zindh, zinda, zindb, ztmu 
     91      REAL(wp) ::   zsto, zjulian, zout,   &  ! temporary scalars 
     92         &          zindh, zinda, zindb, ztmu 
    8893      REAL(wp), DIMENSION(1)                ::   zdept 
    8994      REAL(wp), DIMENSION(jpi,jpj)          ::   zfield 
     
    127132      DO jj = 2 , jpjm1 
    128133         DO ji = 1 , jpim1   ! NO vector opt. 
    129             zindh  = MAX( rzero , SIGN( rone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) ) 
    130             zinda  = MAX( rzero , SIGN( rone , ( 1.0 - frld(ji,jj) ) - 0.10 ) ) 
     134            zindh  = MAX( zzero , SIGN( zone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) ) 
     135            zinda  = MAX( zzero , SIGN( zone , ( 1.0 - frld(ji,jj) ) - 0.10 ) ) 
    131136            zindb  = zindh * zinda 
    132             ztmu   = MAX( 0.5 * rone , ( tmu(ji,jj) + tmu(ji+1,jj) + tmu(ji,jj+1) + tmu(ji+1,jj+1) ) )  
     137            ztmu   = MAX( 0.5 * zone , ( tmu(ji,jj) + tmu(ji+1,jj) + tmu(ji,jj+1) + tmu(ji+1,jj+1) ) )  
    133138            zcmo(ji,jj,1)  = hsnif (ji,jj) 
    134139            zcmo(ji,jj,2)  = hicif (ji,jj) 
     
    138143            zcmo(ji,jj,6)  = fbif  (ji,jj) 
    139144            zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
    140                &                      + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
    141                &                  / ztmu  
     145                                      + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
     146                                  / ztmu  
    142147 
    143148            zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   & 
    144                &                      + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
    145                &                  / ztmu 
     149                                      + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 
     150                                  / ztmu 
    146151            zcmo(ji,jj,9)  = sst_m(ji,jj) 
    147152            zcmo(ji,jj,10) = sss_m(ji,jj) 
     
    195200      !! ** input   :   Namelist namicewri 
    196201      !!------------------------------------------------------------------- 
    197       INTEGER ::   jf      ! dummy loop indices 
     202      INTEGER ::   nf      ! ??? 
    198203      TYPE FIELD  
    199204         CHARACTER(len = 35) :: ztitle  
     
    204209         REAL                :: zcadd         
    205210      END TYPE FIELD 
    206       TYPE(FIELD) ::   field_1 , field_2 , field_3 , field_4 , field_5 , field_6   
    207       TYPE(FIELD) ::   field_7 , field_8 , field_9 , field_10, field_11, field_12 
    208       TYPE(FIELD) ::   field_13, field_14, field_15, field_16, field_17, field_18, field_19 
     211      TYPE(FIELD) ::  & 
     212         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   & 
     213         field_7 , field_8 , field_9 , field_10, field_11, field_12,   & 
     214         field_13, field_14, field_15, field_16, field_17, field_18,   & 
     215         field_19 
    209216      TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield 
    210       !! 
    211       NAMELIST/namiceout/ noumef,                                          & 
    212          field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,       & 
    213          field_7 , field_8 , field_9 , field_10, field_11, field_12,       & 
    214          field_13, field_14, field_15, field_16, field_17, field_18, field_19 
     217 
     218      NAMELIST/namiceout/ noumef, & 
     219         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   & 
     220         field_7 , field_8 , field_9 , field_10, field_11, field_12,   & 
     221         field_13, field_14, field_15, field_16, field_17, field_18,   & 
     222         field_19 
    215223      !!------------------------------------------------------------------- 
    216224 
     
    238246      zfield(19) = field_19 
    239247       
    240       DO jf = 1, noumef 
    241          titn  (jf) = zfield(jf)%ztitle 
    242          nam   (jf) = zfield(jf)%zname 
    243          uni   (jf) = zfield(jf)%zunit 
    244          nc    (jf) = zfield(jf)%znc 
    245          cmulti(jf) = zfield(jf)%zcmulti 
    246          cadd  (jf) = zfield(jf)%zcadd 
     248      DO nf = 1, noumef 
     249         titn  (nf) = zfield(nf)%ztitle 
     250         nam   (nf) = zfield(nf)%zname 
     251         uni   (nf) = zfield(nf)%zunit 
     252         nc    (nf) = zfield(nf)%znc 
     253         cmulti(nf) = zfield(nf)%zcmulti 
     254         cadd  (nf) = zfield(nf)%zcadd 
    247255      END DO 
    248256 
     
    254262         WRITE(numout,*) '           title                            name     unit      Saving (1/0) ',   & 
    255263            &            '    multiplicative constant       additive constant ' 
    256          DO jf = 1 , noumef          
    257             WRITE(numout,*) '   ', titn(jf), '   ', nam(jf),'      ', uni(jf),'  ', nc(jf),'        ', cmulti(jf),   & 
    258                &       '        ', cadd(jf) 
     264         DO nf = 1 , noumef          
     265            WRITE(numout,*) '   ', titn(nf), '   ', nam(nf),'      ', uni(nf),'  ', nc(nf),'        ', cmulti(nf),   & 
     266               &       '        ', cadd(nf) 
    259267         END DO 
    260268      ENDIF 
     
    273281      !!        Used to find errors in the initial state or save the last 
    274282      !!      ocean state in case of abnormal end of a simulation 
     283      !! 
     284      !! History : 
     285      !!   2.0  !  2009-06  (B. Lemaire) 
    275286      !!---------------------------------------------------------------------- 
    276287      INTEGER, INTENT( in ) ::   kt               ! ocean time-step index) 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_3/limtab.F90

    r1855 r1857  
    22   !!====================================================================== 
    33   !!                       ***  MODULE limtab   *** 
    4    !!   LIM-3 ice model : transform 1D (2D) array to a 2D (1D) array 
     4   !!              transform 1D (2D) array to a 2D (1D) table 
    55   !!====================================================================== 
    66#if defined key_lim3 
    77   !!---------------------------------------------------------------------- 
    8    !!   'key_lim3' :                                  LIM 2.0 sea-ice model 
     8   !!   'key_lim3'                                      LIM3 sea-ice model 
    99   !!---------------------------------------------------------------------- 
    1010   !!   tab_2d_1d  : 2-D to 1-D 
    1111   !!   tab_1d_2d  : 1-D to 2-D 
    1212   !!---------------------------------------------------------------------- 
     13   !! * Modules used 
    1314   USE par_kind 
    1415 
     
    1617   PRIVATE 
    1718 
    18    PUBLIC   tab_2d_1d   ! called by lim_thd 
    19    PUBLIC   tab_1d_2d   ! called by lim_thd 
     19   !! * Routine accessibility 
     20   PUBLIC tab_2d_1d  ! called by lim_ther 
     21   PUBLIC tab_1d_2d  ! called by lim_ther 
    2022 
    2123   !!---------------------------------------------------------------------- 
    22    !! NEMO/LIM 3.3,  UCL-LOCEAN-IPSL (2010)  
     24   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)  
    2325   !! $Id$ 
    2426   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    2628CONTAINS 
    2729 
    28    SUBROUTINE tab_2d_1d( kdim1d, ptab1d, ptab2d, kdim2d_i, kdim2d_j, ptab_ind ) 
    29       !!------------------------------------------------------------------- 
    30       INTEGER , INTENT(in   )                                ::   kdim1d 
    31       INTEGER , INTENT(in   )                                ::   kdim2d_i, kdim2d_j 
    32       REAL(wp), INTENT(in   ), DIMENSION(kdim2d_i, kdim2d_j) ::   ptab2d 
    33       INTEGER , INTENT(in   ), DIMENSION(kdim1d)             ::   ptab_ind 
    34       REAL(wp), INTENT(  out), DIMENSION(kdim1d)             ::   ptab1d 
    35       !! 
    36       INTEGER ::   jn , jid, jjd   ! dummy loop indices 
    37       !!------------------------------------------------------------------- 
    38       !  
    39       DO jn = 1, kdim1d 
    40          jid = MOD( ptab_ind(jn) - 1, kdim2d_i ) + 1 
    41          jjd = ( ptab_ind(jn) - 1 ) / kdim2d_i + 1 
    42          ptab1d(jn) = ptab2d(jid,jjd) 
    43       END DO  
    44       ! 
     30   SUBROUTINE tab_2d_1d ( ndim1d, tab1d, tab2d, ndim2d_x, ndim2d_y, tab_ind ) 
     31 
     32      INTEGER, INTENT(in) :: & 
     33         ndim1d, ndim2d_x, ndim2d_y 
     34 
     35      REAL(wp), DIMENSION (ndim2d_x, ndim2d_y), INTENT(in) ::  & 
     36         tab2d 
     37 
     38      INTEGER, DIMENSION ( ndim1d), INTENT ( in) :: & 
     39         tab_ind 
     40 
     41      REAL(wp), DIMENSION(ndim1d), INTENT ( out) ::  &  
     42         tab1d 
     43 
     44      INTEGER ::  & 
     45         jn , jid, jjd 
     46 
     47      DO jn = 1, ndim1d 
     48         jid        = MOD( tab_ind(jn) - 1, ndim2d_x ) + 1 
     49         jjd        = ( tab_ind(jn) - 1 ) / ndim2d_x + 1 
     50         tab1d( jn) = tab2d( jid, jjd) 
     51      END DO 
     52 
    4553   END SUBROUTINE tab_2d_1d 
    4654 
    4755 
    48    SUBROUTINE tab_1d_2d( kdim1d, ptab2d, ptab_ind, ptab1d, kdim2d_i, kdim2d_j ) 
    49       !!------------------------------------------------------------------- 
    50       INTEGER , INTENT(in   )                                ::   kdim1d 
    51       INTEGER , INTENT(in   )                                ::   kdim2d_i, kdim2d_j 
    52       INTEGER , INTENT(in   ), DIMENSION(kdim1d)             ::   ptab_ind 
    53       REAL(wp), INTENT(in   ), DIMENSION(kdim1d)             ::   ptab1d 
    54       REAL(wp), INTENT(  out), DIMENSION(kdim2d_i, kdim2d_j) ::   ptab2d 
    55       !! 
    56       INTEGER ::   jn, jid, jjd   ! dummy loop indices 
    57       !!------------------------------------------------------------------- 
    58       ! 
    59       DO jn = 1, kdim1d 
    60          jid = MOD( ptab_ind(jn) - 1, kdim2d_i) + 1 
    61          jjd =    ( ptab_ind(jn) - 1 ) / kdim2d_i  + 1 
    62          ptab2d(jid,jjd) = ptab1d(jn) 
     56   SUBROUTINE tab_1d_2d ( ndim1d, tab2d, tab_ind, tab1d, ndim2d_x, ndim2d_y ) 
     57 
     58      INTEGER, INTENT ( in) :: & 
     59         ndim1d, ndim2d_x, ndim2d_y 
     60 
     61      INTEGER, DIMENSION (ndim1d) , INTENT (in) :: & 
     62         tab_ind 
     63 
     64      REAL(wp), DIMENSION(ndim1d), INTENT (in) ::  & 
     65         tab1d   
     66 
     67      REAL(wp), DIMENSION (ndim2d_x, ndim2d_y), INTENT ( out) :: & 
     68         tab2d 
     69 
     70      INTEGER :: & 
     71         jn, jid, jjd 
     72 
     73      DO jn = 1, ndim1d 
     74         jid             = MOD( tab_ind(jn) - 1, ndim2d_x) + 1 
     75         jjd             =    ( tab_ind(jn) - 1 ) / ndim2d_x  + 1 
     76         tab2d(jid, jjd) = tab1d( jn) 
    6377      END DO 
    64       ! 
     78 
    6579   END SUBROUTINE tab_1d_2d 
    6680 
    67 #else 
    68    !!---------------------------------------------------------------------- 
    69    !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    70    !!---------------------------------------------------------------------- 
    7181#endif 
    72  
    73    !!====================================================================== 
    7482END MODULE limtab 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r1855 r1857  
    4141   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   alb_ice   !: albedo of ice 
    4242 
     43   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau_ice    !: u-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2] 
     44   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau_ice    !: v-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2] 
    4345   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr1_i0      !: 1st fraction of sol. rad. which penetrate inside the ice cover 
    4446   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr2_i0      !: 2nd fraction of sol. rad. which penetrate inside the ice cover 
    4547   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_ice     !: solid freshwater budget over ice: sublivation - snow 
    46    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau_ice    !: u-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2] 
    47    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau_ice    !: v-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2] 
    4848 
    4949# if defined key_lim3 
Note: See TracChangeset for help on using the changeset viewer.