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 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 – NEMO

Ignore:
Timestamp:
2015-12-01T16:35:30+01:00 (8 years ago)
Author:
timgraham
Message:

Upgraded branch to r5518 of trunk (v3.6 stable revision)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r4624 r5965  
    66   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code 
    77   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
    8    !!             -   ! 2012    (C. Rousset) add par_oce (for jp_sal)...bug? 
     8   !!             -   ! 2014    (C. Rousset) add N/S initializations 
    99   !!---------------------------------------------------------------------- 
    1010#if defined key_lim3 
     
    2222   USE eosbn2           ! equation of state 
    2323   USE ice              ! sea-ice variables 
    24    USE par_ice          ! ice parameters 
    2524   USE par_oce          ! ocean parameters 
    2625   USE dom_ice          ! sea-ice domain 
    2726   USE in_out_manager   ! I/O manager 
    28    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    2927   USE lib_mpp          ! MPP library 
    3028   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     
    3634   PUBLIC   lim_istate      ! routine called by lim_init.F90 
    3735 
    38    !! * Module variables 
    3936   !                          !!** init namelist (namiceini) ** 
    40    REAL(wp) ::   ttest   ! threshold water temperature for initial sea ice 
    41    REAL(wp) ::   hninn   ! initial snow thickness in the north 
    42    REAL(wp) ::   hnins   ! initial snow thickness in the south 
    43    REAL(wp) ::   hginn   ! initial ice thickness in the north 
    44    REAL(wp) ::   hgins   ! initial ice thickness in the south 
    45    REAL(wp) ::   aginn   ! initial leads area in the north 
    46    REAL(wp) ::   agins   ! initial leads area in the south 
    47    REAL(wp) ::   sinn    ! initial salinity  
    48    REAL(wp) ::   sins   
    49  
     37   REAL(wp) ::   rn_thres_sst   ! threshold water temperature for initial sea ice 
     38   REAL(wp) ::   rn_hts_ini_n   ! initial snow thickness in the north 
     39   REAL(wp) ::   rn_hts_ini_s   ! initial snow thickness in the south 
     40   REAL(wp) ::   rn_hti_ini_n   ! initial ice thickness in the north 
     41   REAL(wp) ::   rn_hti_ini_s   ! initial ice thickness in the south 
     42   REAL(wp) ::   rn_ati_ini_n   ! initial leads area in the north 
     43   REAL(wp) ::   rn_ati_ini_s   ! initial leads area in the south 
     44   REAL(wp) ::   rn_smi_ini_n   ! initial salinity  
     45   REAL(wp) ::   rn_smi_ini_s   ! initial salinity 
     46   REAL(wp) ::   rn_tmi_ini_n   ! initial temperature 
     47   REAL(wp) ::   rn_tmi_ini_s   ! initial temperature 
     48 
     49   LOGICAL  ::  ln_iceini    ! initialization or not 
    5050   !!---------------------------------------------------------------------- 
    5151   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008) 
     
    5353   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5454   !!---------------------------------------------------------------------- 
    55  
    5655CONTAINS 
    5756 
     
    7776      !! 
    7877      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even 
    79       !!              where there is no ice (clem: I do not know why but it is mandatory)  
     78      !!              where there is no ice (clem: I do not know why, is it mandatory?)  
    8079      !! 
    8180      !! History : 
     
    8786      !! * Local variables 
    8887      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices 
    89       REAL(wp)   :: epsi20, ztmelts, zdh 
     88      REAL(wp)   :: ztmelts, zdh 
    9089      INTEGER    :: i_hemis, i_fill, jl0   
    9190      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
    92       REAL(wp), POINTER, DIMENSION(:)     :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini 
    93       REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, za_i_ini, zv_i_ini 
    94       REAL(wp), POINTER, DIMENSION(:,:)   :: zidto    ! ice indicator 
     91      REAL(wp), POINTER, DIMENSION(:)     :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini 
     92      REAL(wp), POINTER, DIMENSION(:,:)   :: zh_i_ini, za_i_ini, zv_i_ini 
     93      REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator 
    9594      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index 
    9695      !-------------------------------------------------------------------- 
    9796 
    98       CALL wrk_alloc( jpi, jpj, zidto ) 
     97      CALL wrk_alloc( jpi, jpj, zswitch ) 
    9998      CALL wrk_alloc( jpi, jpj, zhemis ) 
    100       CALL wrk_alloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini ) 
    101       CALL wrk_alloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 
    102  
    103       epsi20   = 1.0e-20 
     99      CALL wrk_alloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini ) 
     100      CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 
     101 
    104102      IF(lwp) WRITE(numout,*) 
    105103      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 
     
    112110      CALL lim_istate_init     !  reading the initials parameters of the ice 
    113111 
    114 !!gm  in lim2  the initialisation if only done if required in the namelist : 
    115 !!gm      IF( .NOT. ln_limini ) THEN 
    116 !!gm  this should be added in lim3 namelist... 
     112      ! surface temperature 
     113      DO jl = 1, jpl ! loop over categories 
     114         t_su  (:,:,jl) = rt0 * tmask(:,:,1) 
     115         tn_ice(:,:,jl) = rt0 * tmask(:,:,1) 
     116      END DO 
     117 
     118      ! basal temperature (considered at freezing point) 
     119      t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1)  
     120 
     121      IF( ln_iceini ) THEN 
    117122 
    118123      !-------------------------------------------------------------------- 
    119124      ! 2) Basal temperature, ice mask and hemispheric index 
    120125      !-------------------------------------------------------------------- 
    121  
    122       ! Basal temperature is set to the freezing point of seawater in Celsius 
    123       t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
    124126 
    125127      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
    126128         DO ji = 1, jpi 
    127             IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0._wp      ! no ice 
    128             ELSE                                                     ;   zidto(ji,jj) = 1._wp      !    ice 
     129            IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN  
     130               zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice 
     131            ELSE                                                                                    
     132               zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice 
    129133            ENDIF 
    130134         END DO 
    131135      END DO 
    132136 
    133       t_bo(:,:) = t_bo(:,:) + rt0                          ! conversion to Kelvin 
    134137 
    135138      ! Hemispheric index 
    136       ! MV 2011 new initialization 
    137139      DO jj = 1, jpj 
    138140         DO ji = 1, jpi 
     
    144146         END DO 
    145147      END DO 
    146       ! END MV 2011 new initialization 
    147148 
    148149      !-------------------------------------------------------------------- 
     
    153154      ! 3.1) Hemisphere-dependent arrays 
    154155      !----------------------------- 
    155       ! assign initial thickness, concentration, snow depth and salinity to 
    156       ! an hemisphere-dependent array 
    157       zhm_i_ini(1) = hginn ; zhm_i_ini(2) = hgins  ! ice thickness 
    158       zat_i_ini(1) = aginn ; zat_i_ini(2) = agins  ! ice concentration 
    159       zvt_i_ini(:) = zhm_i_ini(:) * zat_i_ini(:)   ! ice volume 
    160       zhm_s_ini(1) = hninn ; zhm_s_ini(2) = hnins  ! snow depth 
    161       zsm_i_ini(1) = sinn  ; zsm_i_ini(2) = sins   ! bulk ice salinity 
     156      ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 
     157      zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s  ! ice thickness 
     158      zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s  ! snow depth 
     159      zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s  ! ice concentration 
     160      zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s  ! bulk ice salinity 
     161      ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s  ! temperature (ice and snow) 
     162 
     163      zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:)   ! ice volume 
    162164 
    163165      !--------------------------------------------------------------------- 
     
    183185            ! *** 1 category to fill 
    184186            IF ( i_fill .EQ. 1 ) THEN 
    185                zht_i_ini(1,i_hemis)       = zhm_i_ini(i_hemis) 
    186                za_i_ini(1,i_hemis)        = zat_i_ini(i_hemis) 
    187                zht_i_ini(2:jpl,i_hemis)   = 0._wp 
    188                za_i_ini(2:jpl,i_hemis)    = 0._wp 
     187               zh_i_ini(1,i_hemis)       = zht_i_ini(i_hemis) 
     188               za_i_ini(1,i_hemis)       = zat_i_ini(i_hemis) 
     189               zh_i_ini(2:jpl,i_hemis)   = 0._wp 
     190               za_i_ini(2:jpl,i_hemis)   = 0._wp 
    189191            ELSE 
    190192 
    191             ! *** >1 categores to fill 
    192             !--- Ice thicknesses in the i_fill - 1 first categories 
     193               ! *** >1 categores to fill 
     194               !--- Ice thicknesses in the i_fill - 1 first categories 
    193195               DO jl = 1, i_fill - 1 
    194                   zht_i_ini(jl,i_hemis)    = 0.5 * ( hi_max(jl) + hi_max(jl-1) ) 
    195                END DO 
    196  
    197             !--- jl0: most likely index where cc will be maximum 
     196                  zh_i_ini(jl,i_hemis) = hi_mean(jl) 
     197               END DO 
     198                
     199               !--- jl0: most likely index where cc will be maximum 
    198200               DO jl = 1, jpl 
    199                   IF ( ( zhm_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. & 
    200                        ( zhm_i_ini(i_hemis) .LE. hi_max(jl)   ) ) THEN 
     201                  IF ( ( zht_i_ini(i_hemis) > hi_max(jl-1) ) .AND. & 
     202                     & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN 
    201203                     jl0 = jl 
    202204                  ENDIF 
    203205               END DO 
    204206               jl0 = MIN(jl0, i_fill) 
    205  
    206             !--- Concentrations 
     207                
     208               !--- Concentrations 
    207209               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) 
    208210               DO jl = 1, i_fill - 1 
    209211                  IF ( jl .NE. jl0 ) THEN 
    210                      zsigma               = 0.5 * zhm_i_ini(i_hemis) 
    211                      zarg                 = ( zht_i_ini(jl,i_hemis) - zhm_i_ini(i_hemis) ) / zsigma 
     212                     zsigma               = 0.5 * zht_i_ini(i_hemis) 
     213                     zarg                 = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma 
    212214                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) 
    213215                  ENDIF 
    214                END DO  
    215  
     216               END DO 
     217                
    216218               zA = 0. ! sum of the areas in the jpl categories  
    217219               DO jl = 1, i_fill - 1 
     
    221223               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
    222224          
    223             !--- Ice thickness in the last category 
     225               !--- Ice thickness in the last category 
    224226               zV = 0. ! sum of the volumes of the N-1 categories 
    225227               DO jl = 1, i_fill - 1 
    226                   zV = zV + za_i_ini(jl,i_hemis)*zht_i_ini(jl,i_hemis) 
    227                END DO 
    228                zht_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis)  
    229                IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
    230  
    231             !--- volumes 
    232                zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zht_i_ini(:,i_hemis) 
     228                  zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis) 
     229               END DO 
     230               zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis)  
     231               IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
     232 
     233               !--- volumes 
     234               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis) 
    233235               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
    234236 
     
    262264 
    263265            ! Test 3: thickness of the last category is in-bounds ? 
    264             IF ( zht_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN 
     266            IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN 
    265267               ztest_3 = 1 
    266268            ELSE 
    267269               ! this write is useful 
    268                IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,i_hemis) = ', & 
    269                zht_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 
     270               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', & 
     271               zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 
    270272               ztest_3 = 0 
    271273            ENDIF 
     
    288290 
    289291      IF(lwp) THEN  
    290          WRITE(numout,*), ' ztests : ', ztests 
     292         WRITE(numout,*) ' ztests : ', ztests 
    291293         IF ( ztests .NE. 4 ) THEN 
    292294            WRITE(numout,*) 
    293             WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
    294             WRITE(numout,*), ' !!!! RED ALERT                  !!! ' 
    295             WRITE(numout,*), ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!' 
    296             WRITE(numout,*), ' !!!! Something is wrong in the LIM3 initialization procedure ' 
    297             WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     295            WRITE(numout,*) ' !!!! ALERT                  !!! ' 
     296            WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 
    298297            WRITE(numout,*) 
    299             WRITE(numout,*), ' *** ztests is not equal to 4 ' 
    300             WRITE(numout,*), ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
    301             WRITE(numout,*), ' zat_i_ini : ', zat_i_ini(i_hemis) 
    302             WRITE(numout,*), ' zhm_i_ini : ', zhm_i_ini(i_hemis) 
     298            WRITE(numout,*) ' *** ztests is not equal to 4 ' 
     299            WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
     300            WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis) 
     301            WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis) 
    303302         ENDIF ! ztests .NE. 4 
    304303      ENDIF 
     
    314313         DO jj = 1, jpj 
    315314            DO ji = 1, jpi 
    316                a_i(ji,jj,jl)   = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
    317                ht_i(ji,jj,jl)  = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj))  ! ice thickness 
    318                ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zhm_s_ini( zhemis(ji,jj) ) / zhm_i_ini( zhemis(ji,jj) ) )  ! snow depth 
    319                sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min ! salinity 
    320                o_i(ji,jj,jl)   = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age 
    321                t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * 270.0 ! surf temp 
     315               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
     316               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness 
     317               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth 
     318               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))    ! salinity 
     319               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day) 
     320               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 
    322321 
    323322               ! This case below should not be used if (ht_s/ht_i) is ok in namelist 
     
    327326               ! recompute ht_i, ht_s avoiding out of bounds values 
    328327               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 
    329                ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) 
     328               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 
    330329 
    331330               ! ice volume, salt content, age content 
     
    334333               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
    335334               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content 
    336             END DO ! ji 
    337          END DO ! jj 
    338       END DO ! jl 
     335            END DO 
     336         END DO 
     337      END DO 
    339338 
    340339      ! Snow temperature and heat content 
     
    343342            DO jj = 1, jpj 
    344343               DO ji = 1, jpi 
    345                    t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * rtt 
     344                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 
    346345                   ! Snow energy of melting 
    347                    e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    348                    ! Change dimensions 
    349                    e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    350                    ! Multiply by volume, so that heat content in 10^9 Joules 
    351                    e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 
    352                END DO ! ji 
    353             END DO ! jj 
    354          END DO ! jl 
    355       END DO ! jk 
     346                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
     347 
     348                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
     349                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
     350               END DO 
     351            END DO 
     352         END DO 
     353      END DO 
    356354 
    357355      ! Ice salinity, temperature and heat content 
     
    360358            DO jj = 1, jpj 
    361359               DO ji = 1, jpi 
    362                    t_i(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1._wp - zidto(ji,jj) ) * rtt  
    363                    s_i(ji,jj,jk,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min 
    364                    ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
     360                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0  
     361                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin 
     362                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 
    365363 
    366364                   ! heat content per unit volume 
    367                    e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    368                       +   lfus    * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 
    369                       -   rcp     * ( ztmelts - rtt ) ) 
    370  
    371                    ! Correct dimensions to avoid big values 
    372                    e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    373  
    374                    ! Mutliply by ice volume, and divide by number of layers  
    375                    ! to get heat content in 10^9 J 
    376                    e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i 
    377                END DO ! ji 
    378             END DO ! jj 
    379          END DO ! jl 
    380       END DO ! jk 
    381  
     365                   e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     366                      +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
     367                      -   rcp     * ( ztmelts - rt0 ) ) 
     368 
     369                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
     370                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
     371               END DO 
     372            END DO 
     373         END DO 
     374      END DO 
     375 
     376      tn_ice (:,:,:) = t_su (:,:,:) 
     377 
     378      ELSE  
     379         ! if ln_iceini=false 
     380         a_i  (:,:,:) = 0._wp 
     381         v_i  (:,:,:) = 0._wp 
     382         v_s  (:,:,:) = 0._wp 
     383         smv_i(:,:,:) = 0._wp 
     384         oa_i (:,:,:) = 0._wp 
     385         ht_i (:,:,:) = 0._wp 
     386         ht_s (:,:,:) = 0._wp 
     387         sm_i (:,:,:) = 0._wp 
     388         o_i  (:,:,:) = 0._wp 
     389 
     390         e_i(:,:,:,:) = 0._wp 
     391         e_s(:,:,:,:) = 0._wp 
     392 
     393         DO jl = 1, jpl 
     394            DO jk = 1, nlay_i 
     395               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 
     396            END DO 
     397            DO jk = 1, nlay_s 
     398               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 
     399            END DO 
     400         END DO 
     401       
     402      ENDIF ! ln_iceini 
     403       
     404      at_i (:,:) = 0.0_wp 
     405      DO jl = 1, jpl 
     406         at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 
     407      END DO 
     408      ! 
    382409      !-------------------------------------------------------------------- 
    383410      ! 4) Global ice variables for output diagnostics                    |  
    384411      !-------------------------------------------------------------------- 
    385       fsbbq (:,:)     = 0._wp 
    386412      u_ice (:,:)     = 0._wp 
    387413      v_ice (:,:)     = 0._wp 
     
    390416      stress12_i(:,:) = 0._wp 
    391417 
    392 # if defined key_coupled 
    393       albege(:,:)   = 0.8 * tms(:,:) 
    394 # endif 
    395  
    396418      !-------------------------------------------------------------------- 
    397419      ! 5) Moments for advection 
     
    428450      sxyage (:,:,:)  = 0._wp 
    429451 
    430       !-------------------------------------------------------------------- 
    431       ! 6) Lateral boundary conditions                                    |  
    432       !-------------------------------------------------------------------- 
    433  
    434       DO jl = 1, jpl 
    435  
    436          CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. ) 
    437          CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. ) 
    438          CALL lbc_lnk( v_s(:,:,jl)  , 'T', 1. ) 
    439          CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. ) 
    440          CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. ) 
    441  
    442          CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. ) 
    443          CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. ) 
    444          CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. ) 
    445          CALL lbc_lnk( o_i(:,:,jl)  , 'T', 1. ) 
    446          CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. ) 
    447          DO jk = 1, nlay_s 
    448             CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. ) 
    449             CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. ) 
    450          END DO 
    451          DO jk = 1, nlay_i 
    452             CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. ) 
    453             CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. ) 
    454          END DO 
    455          ! 
    456          a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl) 
    457       END DO 
    458        
    459       at_i (:,:) = 0.0_wp 
    460       DO jl = 1, jpl 
    461          at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 
    462       END DO 
    463  
    464       CALL lbc_lnk( at_i , 'T', 1. ) 
    465       at_i(:,:) = tms(:,:) * at_i(:,:)                       ! put 0 over land 
    466       ! 
    467       CALL lbc_lnk( fsbbq  , 'T', 1. ) 
    468       ! 
    469       !-------------------------------------------------------------------- 
    470       ! 6) ????                                                           |  
    471       !-------------------------------------------------------------------- 
    472       tn_ice (:,:,:) = t_su (:,:,:) 
    473  
    474       CALL wrk_dealloc( jpi, jpj, zidto ) 
     452 
     453      CALL wrk_dealloc( jpi, jpj, zswitch ) 
    475454      CALL wrk_dealloc( jpi, jpj, zhemis ) 
    476       CALL wrk_dealloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini ) 
    477       CALL wrk_dealloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 
     455      CALL wrk_dealloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini ) 
     456      CALL wrk_dealloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 
    478457 
    479458   END SUBROUTINE lim_istate 
     
    495474      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization 
    496475      !!----------------------------------------------------------------------------- 
    497       NAMELIST/namiceini/ ttest, hninn, hnins, hginn, hgins, aginn, agins, sinn, sins 
    498       ! 
     476      NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s,  & 
     477         &                                      rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s 
    499478      INTEGER :: ios                 ! Local integer output status for namelist read 
    500479      !!----------------------------------------------------------------------------- 
     
    516495         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 
    517496         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    518          WRITE(numout,*) '   threshold water temp. for initial sea-ice    ttest      = ', ttest 
    519          WRITE(numout,*) '   initial snow thickness in the north          hninn      = ', hninn 
    520          WRITE(numout,*) '   initial snow thickness in the south          hnins      = ', hnins  
    521          WRITE(numout,*) '   initial ice thickness  in the north          hginn      = ', hginn 
    522          WRITE(numout,*) '   initial ice thickness  in the south          hgins      = ', hgins 
    523          WRITE(numout,*) '   initial ice concentr.  in the north          aginn      = ', aginn 
    524          WRITE(numout,*) '   initial ice concentr.  in the north          agins      = ', agins 
    525          WRITE(numout,*) '   initial  ice salinity  in the north          sinn       = ', sinn 
    526          WRITE(numout,*) '   initial  ice salinity  in the south          sins       = ', sins 
     497         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini 
     498         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst 
     499         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n 
     500         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s  
     501         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n 
     502         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s 
     503         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n 
     504         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s 
     505         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n 
     506         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s 
     507         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n 
     508         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s 
    527509      ENDIF 
    528510 
Note: See TracChangeset for help on using the changeset viewer.