MODULE limistate !!====================================================================== !! *** MODULE limistate *** !! Initialisation of diagnostics ice variables !!====================================================================== !! History : 2.0 ! 2004-01 (C. Ethe, G. Madec) Original code !! 4.0 ! 2011-02 (G. Madec) dynamical allocation !! - ! 2012 (C. Rousset) add par_oce (for jp_sal)...bug? !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' : LIM3 sea-ice model !!---------------------------------------------------------------------- !! lim_istate : Initialisation of diagnostics ice variables !! lim_istate_init : initialization of ice state and namelist read !!---------------------------------------------------------------------- USE phycst ! physical constant USE oce ! dynamics and tracers variables USE dom_oce ! ocean domain USE sbc_oce ! Surface boundary condition: ocean fields USE sbc_ice ! Surface boundary condition: ice fields USE eosbn2 ! equation of state USE ice ! sea-ice variables USE par_ice ! ice parameters USE par_oce ! ocean parameters USE dom_ice ! sea-ice domain USE in_out_manager ! I/O manager USE lbclnk ! lateral boundary condition - MPP exchanges USE lib_mpp ! MPP library USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) USE wrk_nemo ! work arrays IMPLICIT NONE PRIVATE PUBLIC lim_istate ! routine called by lim_init.F90 !! * Module variables REAL(wp) :: & !!! ** init namelist (namiceini) ** ttest = 2.0 , & ! threshold water temperature for initial sea ice hninn = 0.5 , & ! initial snow thickness in the north hnins = 0.1 , & ! initial snow thickness in the south hginn = 2.5 , & ! initial ice thickness in the north hgins = 1.0 , & ! initial ice thickness in the south aginn = 0.7 , & ! initial leads area in the north agins = 0.7 , & ! initial leads area in the south sinn = 6.301 , & ! initial salinity sins = 6.301 !!---------------------------------------------------------------------- !! LIM 3.0, UCL-LOCEAN-IPSL (2008) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_istate !!------------------------------------------------------------------- !! *** ROUTINE lim_istate *** !! !! ** Purpose : defined the sea-ice initial state !! !! ** Method : !! This routine will put some ice where ocean !! is at the freezing point, then fill in ice !! state variables using prescribed initial !! values in the namelist !! !! ** Steps : !! 1) Read namelist !! 2) Basal temperature; ice and hemisphere masks !! 3) Fill in the ice thickness distribution using gaussian !! 4) Fill in space-dependent arrays for state variables !! 5) Diagnostic arrays !! 6) Lateral boundary conditions !! !! ** Notes : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even !! where there is no ice (clem: I do not know why but it is mandatory) !! !! History : !! 2.0 ! 01-04 (C. Ethe, G. Madec) Original code !! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats !! 4.0 ! 09-11 (M. Vancoppenolle) Enhanced version for ice cats !!-------------------------------------------------------------------- !! * Local variables INTEGER :: ji, jj, jk, jl ! dummy loop indices REAL(wp) :: epsi20, ztmelts, zdh INTEGER :: i_hemis, i_fill, jl0 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv REAL(wp), POINTER, DIMENSION(:) :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, za_i_ini, zv_i_ini REAL(wp), POINTER, DIMENSION(:,:) :: zidto ! ice indicator INTEGER, POINTER, DIMENSION(:,:) :: zhemis ! hemispheric index !-------------------------------------------------------------------- CALL wrk_alloc( jpi, jpj, zidto ) CALL wrk_alloc( jpi, jpj, zhemis ) CALL wrk_alloc( jpl, 2, zht_i_ini, za_i_ini, zv_i_ini ) CALL wrk_alloc( 2, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) epsi20 = 1.0e-20 IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' !-------------------------------------------------------------------- ! 1) Read namelist !-------------------------------------------------------------------- CALL lim_istate_init ! reading the initials parameters of the ice !!gm in lim2 the initialisation if only done if required in the namelist : !!gm IF( .NOT. ln_limini ) THEN !!gm this should be added in lim3 namelist... !-------------------------------------------------------------------- ! 2) Basal temperature, ice mask and hemispheric index !-------------------------------------------------------------------- ! Basal temperature is set to the freezing point of seawater in Celsius t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius] DO jj = 1, jpj ! ice if sst <= t-freez + ttest DO ji = 1, jpi IF( tsn(ji,jj,1,jp_tem) - t_bo(ji,jj) >= ttest ) THEN ; zidto(ji,jj) = 0._wp ! no ice ELSE ; zidto(ji,jj) = 1._wp ! ice ENDIF END DO END DO t_bo(:,:) = t_bo(:,:) + rt0 ! conversion to Kelvin ! Hemispheric index ! MV 2011 new initialization DO jj = 1, jpj DO ji = 1, jpi IF( fcor(ji,jj) >= 0._wp ) THEN zhemis(ji,jj) = 1 ! Northern hemisphere ELSE zhemis(ji,jj) = 2 ! Southern hemisphere ENDIF END DO END DO ! END MV 2011 new initialization !-------------------------------------------------------------------- ! 3) Initialization of sea ice state variables !-------------------------------------------------------------------- !----------------------------- ! 3.1) Hemisphere-dependent arrays !----------------------------- ! assign initial thickness, concentration, snow depth and salinity to ! an hemisphere-dependent array zhm_i_ini(1) = hginn ; zhm_i_ini(2) = hgins ! ice thickness zat_i_ini(1) = aginn ; zat_i_ini(2) = agins ! ice concentration zvt_i_ini(:) = zhm_i_ini(:) * zat_i_ini(:) ! ice volume zhm_s_ini(1) = hninn ; zhm_s_ini(2) = hnins ! snow depth zsm_i_ini(1) = sinn ; zsm_i_ini(2) = sins ! bulk ice salinity !--------------------------------------------------------------------- ! 3.2) Distribute ice concentration and thickness into the categories !--------------------------------------------------------------------- ! a gaussian distribution for ice concentration is used ! then we check whether the distribution fullfills ! volume and area conservation, positivity and ice categories bounds DO i_hemis = 1, 2 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 ! note for the great nemo engineers: ! only very few of the WRITE statements are necessary for the reference version ! they were one day useful, but now i personally doubt of their ! potential for bringing anything useful DO i_fill = jpl, 1, -1 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN !---------------------------- ! fill the i_fill categories !---------------------------- ! *** 1 category to fill IF ( i_fill .EQ. 1 ) THEN zht_i_ini(1,i_hemis) = zhm_i_ini(i_hemis) za_i_ini(1,i_hemis) = zat_i_ini(i_hemis) zht_i_ini(2:jpl,i_hemis) = 0._wp za_i_ini(2:jpl,i_hemis) = 0._wp ELSE ! *** >1 categores to fill !--- Ice thicknesses in the i_fill - 1 first categories DO jl = 1, i_fill - 1 zht_i_ini(jl,i_hemis) = 0.5 * ( hi_max(jl) + hi_max(jl-1) ) END DO !--- jl0: most likely index where cc will be maximum DO jl = 1, jpl IF ( ( zhm_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. & ( zhm_i_ini(i_hemis) .LE. hi_max(jl) ) ) THEN jl0 = jl ENDIF END DO jl0 = MIN(jl0, i_fill) !--- Concentrations za_i_ini(jl0,i_hemis) = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) DO jl = 1, i_fill - 1 IF ( jl .NE. jl0 ) THEN zsigma = 0.5 * zhm_i_ini(i_hemis) zarg = ( zht_i_ini(jl,i_hemis) - zhm_i_ini(i_hemis) ) / zsigma za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) ENDIF END DO zA = 0. ! sum of the areas in the jpl categories DO jl = 1, i_fill - 1 zA = zA + za_i_ini(jl,i_hemis) END DO za_i_ini(i_fill,i_hemis) = zat_i_ini(i_hemis) - zA ! ice conc in the last category IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp !--- Ice thickness in the last category zV = 0. ! sum of the volumes of the N-1 categories DO jl = 1, i_fill - 1 zV = zV + za_i_ini(jl,i_hemis)*zht_i_ini(jl,i_hemis) END DO zht_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, i_hemis) = 0._wp !--- volumes zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zht_i_ini(:,i_hemis) IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp ENDIF ! i_fill !--------------------- ! Compatibility tests !--------------------- ! Test 1: area conservation zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons ) IF ( zconv .LT. 1.0e-6 ) THEN ztest_1 = 1 ELSE ! this write is useful IF(lwp) WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) ztest_1 = 0 ENDIF ! Test 2: volume conservation zV_cons = SUM(zv_i_ini(:,i_hemis)) zconv = ABS(zvt_i_ini(i_hemis) - zV_cons) IF ( zconv .LT. 1.0e-6 ) THEN ztest_2 = 1 ELSE ! this write is useful IF(lwp) WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & ' zvt_i_ini = ', zvt_i_ini(i_hemis) ztest_2 = 0 ENDIF ! Test 3: thickness of the last category is in-bounds ? IF ( zht_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN ztest_3 = 1 ELSE ! this write is useful IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,i_hemis) = ', & zht_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) ztest_3 = 0 ENDIF ! Test 4: positivity of ice concentrations ztest_4 = 1 DO jl = 1, jpl IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN ! this write is useful IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis) ztest_4 = 0 ENDIF END DO ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 END DO ! i_fill IF(lwp) THEN WRITE(numout,*), ' ztests : ', ztests IF ( ztests .NE. 4 ) THEN WRITE(numout,*) WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' WRITE(numout,*), ' !!!! RED ALERT !!! ' WRITE(numout,*), ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!' WRITE(numout,*), ' !!!! Something is wrong in the LIM3 initialization procedure ' WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' WRITE(numout,*) WRITE(numout,*), ' *** ztests is not equal to 4 ' WRITE(numout,*), ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 WRITE(numout,*), ' zat_i_ini : ', zat_i_ini(i_hemis) WRITE(numout,*), ' zhm_i_ini : ', zhm_i_ini(i_hemis) ENDIF ! ztests .NE. 4 ENDIF END DO ! i_hemis !--------------------------------------------------------------------- ! 3.3) Space-dependent arrays for ice state variables !--------------------------------------------------------------------- ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature DO jl = 1, jpl ! loop over categories DO jj = 1, jpj DO ji = 1, jpi a_i(ji,jj,jl) = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj)) ! concentration ht_i(ji,jj,jl) = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj)) ! ice thickness ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zhm_s_ini( zhemis(ji,jj) ) / zhm_i_ini( zhemis(ji,jj) ) ) ! snow depth sm_i(ji,jj,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min ! salinity o_i(ji,jj,jl) = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * 270.0 ! surf temp ! This case below should not be used if (ht_s/ht_i) is ok in namelist ! In case snow load is in excess that would lead to transformation from snow to ice ! Then, transfer the snow excess into the ice (different from limthd_dh) zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) ! recompute ht_i, ht_s avoiding out of bounds values ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) ! ice volume, salt content, age content v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content END DO ! ji END DO ! jj END DO ! jl ! Snow temperature and heat content DO jk = 1, nlay_s DO jl = 1, jpl ! loop over categories DO jj = 1, jpj DO ji = 1, jpi t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * rtt ! Snow energy of melting e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) ! Change dimensions e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac ! Multiply by volume, so that heat content in 10^9 Joules e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s END DO ! ji END DO ! jj END DO ! jl END DO ! jk ! Ice salinity, temperature and heat content DO jk = 1, nlay_i DO jl = 1, jpl ! loop over categories DO jj = 1, jpj DO ji = 1, jpi t_i(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1._wp - zidto(ji,jj) ) * rtt s_i(ji,jj,jk,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K ! heat content per unit volume e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & + lfus * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & - rcp * ( ztmelts - rtt ) ) ! Correct dimensions to avoid big values e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac ! Mutliply by ice volume, and divide by number of layers ! to get heat content in 10^9 J e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i END DO ! ji END DO ! jj END DO ! jl END DO ! jk !-------------------------------------------------------------------- ! 4) Global ice variables for output diagnostics | !-------------------------------------------------------------------- fsbbq (:,:) = 0._wp u_ice (:,:) = 0._wp v_ice (:,:) = 0._wp stress1_i(:,:) = 0._wp stress2_i(:,:) = 0._wp stress12_i(:,:) = 0._wp # if defined key_coupled albege(:,:) = 0.8 * tms(:,:) # endif !-------------------------------------------------------------------- ! 5) Moments for advection !-------------------------------------------------------------------- sxopw (:,:) = 0._wp syopw (:,:) = 0._wp sxxopw(:,:) = 0._wp syyopw(:,:) = 0._wp sxyopw(:,:) = 0._wp sxice (:,:,:) = 0._wp ; sxsn (:,:,:) = 0._wp ; sxa (:,:,:) = 0._wp syice (:,:,:) = 0._wp ; sysn (:,:,:) = 0._wp ; sya (:,:,:) = 0._wp sxxice(:,:,:) = 0._wp ; sxxsn(:,:,:) = 0._wp ; sxxa (:,:,:) = 0._wp syyice(:,:,:) = 0._wp ; syysn(:,:,:) = 0._wp ; syya (:,:,:) = 0._wp sxyice(:,:,:) = 0._wp ; sxysn(:,:,:) = 0._wp ; sxya (:,:,:) = 0._wp sxc0 (:,:,:) = 0._wp ; sxe (:,:,:,:)= 0._wp syc0 (:,:,:) = 0._wp ; sye (:,:,:,:)= 0._wp sxxc0 (:,:,:) = 0._wp ; sxxe (:,:,:,:)= 0._wp syyc0 (:,:,:) = 0._wp ; syye (:,:,:,:)= 0._wp sxyc0 (:,:,:) = 0._wp ; sxye (:,:,:,:)= 0._wp sxsal (:,:,:) = 0._wp sysal (:,:,:) = 0._wp sxxsal (:,:,:) = 0._wp syysal (:,:,:) = 0._wp sxysal (:,:,:) = 0._wp sxage (:,:,:) = 0._wp syage (:,:,:) = 0._wp sxxage (:,:,:) = 0._wp syyage (:,:,:) = 0._wp sxyage (:,:,:) = 0._wp !-------------------------------------------------------------------- ! 6) Lateral boundary conditions | !-------------------------------------------------------------------- DO jl = 1, jpl CALL lbc_lnk( a_i(:,:,jl) , 'T', 1. ) CALL lbc_lnk( v_i(:,:,jl) , 'T', 1. ) CALL lbc_lnk( v_s(:,:,jl) , 'T', 1. ) CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. ) CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. ) CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. ) CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. ) CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. ) CALL lbc_lnk( o_i(:,:,jl) , 'T', 1. ) CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. ) DO jk = 1, nlay_s CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. ) CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. ) END DO DO jk = 1, nlay_i CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. ) CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. ) END DO ! a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl) END DO at_i (:,:) = 0.0_wp DO jl = 1, jpl at_i (:,:) = at_i (:,:) + a_i (:,:,jl) END DO CALL lbc_lnk( at_i , 'T', 1. ) at_i(:,:) = tms(:,:) * at_i(:,:) ! put 0 over land ! CALL lbc_lnk( fsbbq , 'T', 1. ) ! !-------------------------------------------------------------------- ! 6) ???? | !-------------------------------------------------------------------- tn_ice (:,:,:) = t_su (:,:,:) CALL wrk_dealloc( jpi, jpj, zidto ) CALL wrk_dealloc( jpi, jpj, zhemis ) CALL wrk_dealloc( jpl, 2, zht_i_ini, za_i_ini, zv_i_ini ) CALL wrk_dealloc( 2, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) END SUBROUTINE lim_istate SUBROUTINE lim_istate_init !!------------------------------------------------------------------- !! *** ROUTINE lim_istate_init *** !! !! ** Purpose : Definition of initial state of the ice !! !! ** Method : Read the namiceini namelist and check the parameter !! values called at the first timestep (nit000) !! !! ** input : !! Namelist namiceini !! !! history : !! 8.5 ! 03-08 (C. Ethe) original code !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization !!----------------------------------------------------------------------------- NAMELIST/namiceini/ ttest, hninn, hnins, hginn, hgins, aginn, agins, sinn, sins !!----------------------------------------------------------------------------- ! Define the initial parameters ! ------------------------- ! Read Namelist namiceini REWIND ( numnam_ice ) READ ( numnam_ice , namiceini ) IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' WRITE(numout,*) '~~~~~~~~~~~~~~~' WRITE(numout,*) ' threshold water temp. for initial sea-ice ttest = ', ttest WRITE(numout,*) ' initial snow thickness in the north hninn = ', hninn WRITE(numout,*) ' initial snow thickness in the south hnins = ', hnins WRITE(numout,*) ' initial ice thickness in the north hginn = ', hginn WRITE(numout,*) ' initial ice thickness in the south hgins = ', hgins WRITE(numout,*) ' initial ice concentr. in the north aginn = ', aginn WRITE(numout,*) ' initial ice concentr. in the north agins = ', agins WRITE(numout,*) ' initial ice salinity in the north sinn = ', sinn WRITE(numout,*) ' initial ice salinity in the south sins = ', sins ENDIF END SUBROUTINE lim_istate_init #else !!---------------------------------------------------------------------- !! Default option : Empty module NO LIM sea-ice model !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_istate ! Empty routine END SUBROUTINE lim_istate #endif !!====================================================================== END MODULE limistate