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.
limistate.F90 in branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 6469

Last change on this file since 6469 was 6469, checked in by clem, 8 years ago

solve issue #1712

  • Property svn:keywords set to Id
File size: 24.0 KB
RevLine 
[825]1MODULE limistate
2   !!======================================================================
3   !!                     ***  MODULE  limistate  ***
4   !!              Initialisation of diagnostics ice variables
5   !!======================================================================
[2528]6   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
[4161]7   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
[4990]8   !!             -   ! 2014    (C. Rousset) add N/S initializations
[2528]9   !!----------------------------------------------------------------------
[825]10#if defined key_lim3
11   !!----------------------------------------------------------------------
[834]12   !!   'key_lim3' :                                    LIM3 sea-ice model
[825]13   !!----------------------------------------------------------------------
14   !!   lim_istate      :  Initialisation of diagnostics ice variables
15   !!   lim_istate_init :  initialization of ice state and namelist read
16   !!----------------------------------------------------------------------
[2528]17   USE phycst           ! physical constant
18   USE oce              ! dynamics and tracers variables
19   USE dom_oce          ! ocean domain
20   USE sbc_oce          ! Surface boundary condition: ocean fields
[4161]21   USE sbc_ice          ! Surface boundary condition: ice fields
[2528]22   USE eosbn2           ! equation of state
23   USE ice              ! sea-ice variables
[4161]24   USE par_oce          ! ocean parameters
[2528]25   USE dom_ice          ! sea-ice domain
[6469]26   USE limvar           ! lim_var_salprof
[2528]27   USE in_out_manager   ! I/O manager
[2715]28   USE lib_mpp          ! MPP library
[4161]29   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[3294]30   USE wrk_nemo         ! work arrays
[825]31
32   IMPLICIT NONE
33   PRIVATE
34
[2528]35   PUBLIC   lim_istate      ! routine called by lim_init.F90
[825]36
[4166]37   !                          !!** init namelist (namiceini) **
[5123]38   REAL(wp) ::   rn_thres_sst   ! threshold water temperature for initial sea ice
39   REAL(wp) ::   rn_hts_ini_n   ! initial snow thickness in the north
40   REAL(wp) ::   rn_hts_ini_s   ! initial snow thickness in the south
41   REAL(wp) ::   rn_hti_ini_n   ! initial ice thickness in the north
42   REAL(wp) ::   rn_hti_ini_s   ! initial ice thickness in the south
43   REAL(wp) ::   rn_ati_ini_n   ! initial leads area in the north
44   REAL(wp) ::   rn_ati_ini_s   ! initial leads area in the south
45   REAL(wp) ::   rn_smi_ini_n   ! initial salinity
46   REAL(wp) ::   rn_smi_ini_s   ! initial salinity
47   REAL(wp) ::   rn_tmi_ini_n   ! initial temperature
48   REAL(wp) ::   rn_tmi_ini_s   ! initial temperature
[4166]49
[5123]50   LOGICAL  ::  ln_iceini    ! initialization or not
[825]51   !!----------------------------------------------------------------------
[4161]52   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
[1156]53   !! $Id$
[4161]54   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[825]55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE lim_istate
59      !!-------------------------------------------------------------------
60      !!                    ***  ROUTINE lim_istate  ***
61      !!
62      !! ** Purpose :   defined the sea-ice initial state
63      !!
[4161]64      !! ** Method  :   
65      !!                This routine will put some ice where ocean
66      !!                is at the freezing point, then fill in ice
67      !!                state variables using prescribed initial
68      !!                values in the namelist           
69      !!
70      !! ** Steps   :   
71      !!                1) Read namelist
72      !!                2) Basal temperature; ice and hemisphere masks
73      !!                3) Fill in the ice thickness distribution using gaussian
74      !!                4) Fill in space-dependent arrays for state variables
75      !!                5) Diagnostic arrays
76      !!                6) Lateral boundary conditions
77      !!
[4333]78      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
[4990]79      !!              where there is no ice (clem: I do not know why, is it mandatory?)
[4333]80      !!
[4161]81      !! History :
82      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
83      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
84      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
85      !!--------------------------------------------------------------------
86
87      !! * Local variables
88      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
[5123]89      REAL(wp)   :: ztmelts, zdh
[4161]90      INTEGER    :: i_hemis, i_fill, jl0 
91      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv
[4688]92      REAL(wp), POINTER, DIMENSION(:)     :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini
93      REAL(wp), POINTER, DIMENSION(:,:)   :: zh_i_ini, za_i_ini, zv_i_ini
94      REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator
[4161]95      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index
[825]96      !--------------------------------------------------------------------
[921]97
[4688]98      CALL wrk_alloc( jpi, jpj, zswitch )
[4161]99      CALL wrk_alloc( jpi, jpj, zhemis )
[4688]100      CALL wrk_alloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
101      CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
[2715]102
[4161]103      IF(lwp) WRITE(numout,*)
104      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
105      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
106
[825]107      !--------------------------------------------------------------------
[4161]108      ! 1) Read namelist
[825]109      !--------------------------------------------------------------------
110
111      CALL lim_istate_init     !  reading the initials parameters of the ice
112
[4688]113      ! surface temperature
114      DO jl = 1, jpl ! loop over categories
[5123]115         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
116         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
[4688]117      END DO
118
[4990]119      ! basal temperature (considered at freezing point)
[5540]120      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
121      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
[4990]122
[5123]123      IF( ln_iceini ) THEN
[4688]124
[825]125      !--------------------------------------------------------------------
[4161]126      ! 2) Basal temperature, ice mask and hemispheric index
[825]127      !--------------------------------------------------------------------
[4990]128
129      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
[1037]130         DO ji = 1, jpi
[5407]131            IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN
[5123]132               zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
[4765]133            ELSE                                                                                   
[5123]134               zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
[1037]135            ENDIF
136         END DO
137      END DO
138
139
[4161]140      ! Hemispheric index
[825]141      DO jj = 1, jpj
142         DO ji = 1, jpi
[2528]143            IF( fcor(ji,jj) >= 0._wp ) THEN   
[4161]144               zhemis(ji,jj) = 1 ! Northern hemisphere
145            ELSE
146               zhemis(ji,jj) = 2 ! Southern hemisphere
147            ENDIF
148         END DO
149      END DO
[825]150
[4161]151      !--------------------------------------------------------------------
152      ! 3) Initialization of sea ice state variables
153      !--------------------------------------------------------------------
[825]154
[4161]155      !-----------------------------
156      ! 3.1) Hemisphere-dependent arrays
157      !-----------------------------
[4688]158      ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
[5123]159      zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s  ! ice thickness
160      zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s  ! snow depth
161      zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s  ! ice concentration
162      zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s  ! bulk ice salinity
163      ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s  ! temperature (ice and snow)
[825]164
[4688]165      zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:)   ! ice volume
166
[4161]167      !---------------------------------------------------------------------
168      ! 3.2) Distribute ice concentration and thickness into the categories
169      !---------------------------------------------------------------------
170      ! a gaussian distribution for ice concentration is used
171      ! then we check whether the distribution fullfills
172      ! volume and area conservation, positivity and ice categories bounds
173      DO i_hemis = 1, 2
[825]174
[4161]175      ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
[825]176
[4161]177      ! note for the great nemo engineers:
178      ! only very few of the WRITE statements are necessary for the reference version
179      ! they were one day useful, but now i personally doubt of their
180      ! potential for bringing anything useful
[825]181
[4161]182      DO i_fill = jpl, 1, -1
183         IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
184            !----------------------------
185            ! fill the i_fill categories
186            !----------------------------
187            ! *** 1 category to fill
188            IF ( i_fill .EQ. 1 ) THEN
[4688]189               zh_i_ini(1,i_hemis)       = zht_i_ini(i_hemis)
190               za_i_ini(1,i_hemis)       = zat_i_ini(i_hemis)
191               zh_i_ini(2:jpl,i_hemis)   = 0._wp
192               za_i_ini(2:jpl,i_hemis)   = 0._wp
[4161]193            ELSE
[825]194
[4688]195               ! *** >1 categores to fill
196               !--- Ice thicknesses in the i_fill - 1 first categories
[4161]197               DO jl = 1, i_fill - 1
[5123]198                  zh_i_ini(jl,i_hemis) = hi_mean(jl)
[4161]199               END DO
[4688]200               
201               !--- jl0: most likely index where cc will be maximum
[825]202               DO jl = 1, jpl
[5123]203                  IF ( ( zht_i_ini(i_hemis) >  hi_max(jl-1) ) .AND. &
204                     & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN
[4161]205                     jl0 = jl
206                  ENDIF
207               END DO
208               jl0 = MIN(jl0, i_fill)
[4688]209               
210               !--- Concentrations
[4161]211               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl))
212               DO jl = 1, i_fill - 1
213                  IF ( jl .NE. jl0 ) THEN
[4688]214                     zsigma               = 0.5 * zht_i_ini(i_hemis)
215                     zarg                 = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma
[4161]216                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2)
217                  ENDIF
[4688]218               END DO
219               
[4161]220               zA = 0. ! sum of the areas in the jpl categories
221               DO jl = 1, i_fill - 1
222                 zA = zA + za_i_ini(jl,i_hemis)
223               END DO
224               za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category
225               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
226         
[4688]227               !--- Ice thickness in the last category
[4161]228               zV = 0. ! sum of the volumes of the N-1 categories
229               DO jl = 1, i_fill - 1
[4688]230                  zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis)
[4161]231               END DO
[4688]232               zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 
233               IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
[825]234
[4688]235               !--- volumes
236               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis)
[4161]237               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
[921]238
[4161]239            ENDIF ! i_fill
[825]240
[4161]241            !---------------------
242            ! Compatibility tests
243            !---------------------
244            ! Test 1: area conservation
245            zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons )
246            IF ( zconv .LT. 1.0e-6 ) THEN
247               ztest_1 = 1
248            ELSE 
249              ! this write is useful
250              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 
251               ztest_1 = 0
252            ENDIF
[825]253
[4161]254            ! Test 2: volume conservation
255            zV_cons = SUM(zv_i_ini(:,i_hemis))
256            zconv = ABS(zvt_i_ini(i_hemis) - zV_cons)
[825]257
[4161]258            IF ( zconv .LT. 1.0e-6 ) THEN
259               ztest_2 = 1
260            ELSE
261              ! this write is useful
262              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
263                            ' zvt_i_ini = ', zvt_i_ini(i_hemis)
264               ztest_2 = 0
265            ENDIF
[825]266
[4161]267            ! Test 3: thickness of the last category is in-bounds ?
[5123]268            IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN
[4161]269               ztest_3 = 1
270            ELSE
271               ! this write is useful
[4688]272               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', &
273               zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
[4161]274               ztest_3 = 0
275            ENDIF
[825]276
[4161]277            ! Test 4: positivity of ice concentrations
278            ztest_4 = 1
279            DO jl = 1, jpl
280               IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 
281                  ! this write is useful
282                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis)
283                  ztest_4 = 0
[825]284               ENDIF
[4161]285            END DO
[825]286
[4161]287         ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
288 
289         ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
[825]290
[4161]291      END DO ! i_fill
[825]292
[4161]293      IF(lwp) THEN
[4908]294         WRITE(numout,*) ' ztests : ', ztests
[4161]295         IF ( ztests .NE. 4 ) THEN
296            WRITE(numout,*)
[4908]297            WRITE(numout,*) ' !!!! ALERT                  !!! '
298            WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
[4161]299            WRITE(numout,*)
[4908]300            WRITE(numout,*) ' *** ztests is not equal to 4 '
301            WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
302            WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis)
303            WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis)
[4161]304         ENDIF ! ztests .NE. 4
305      ENDIF
306     
307      END DO ! i_hemis
[825]308
[4161]309      !---------------------------------------------------------------------
310      ! 3.3) Space-dependent arrays for ice state variables
311      !---------------------------------------------------------------------
[825]312
[4333]313      ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
[4161]314      DO jl = 1, jpl ! loop over categories
315         DO jj = 1, jpj
316            DO ji = 1, jpi
[4688]317               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
[5202]318               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness
[4688]319               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth
[5202]320               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity
321               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day)
[5123]322               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
[4333]323
324               ! This case below should not be used if (ht_s/ht_i) is ok in namelist
325               ! In case snow load is in excess that would lead to transformation from snow to ice
326               ! Then, transfer the snow excess into the ice (different from limthd_dh)
327               zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
328               ! recompute ht_i, ht_s avoiding out of bounds values
329               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
[5123]330               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
[4333]331
332               ! ice volume, salt content, age content
[4161]333               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
334               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
335               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
336               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
[5202]337            END DO
338         END DO
339      END DO
[825]340
[6469]341      ! for constant salinity in time
342      IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
343         CALL lim_var_salprof
344         smv_i = sm_i * v_i
345      ENDIF
346
[4161]347      ! Snow temperature and heat content
348      DO jk = 1, nlay_s
349         DO jl = 1, jpl ! loop over categories
350            DO jj = 1, jpj
351               DO ji = 1, jpi
[5123]352                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0
[4161]353                   ! Snow energy of melting
[5123]354                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
355
356                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
357                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
[5202]358               END DO
359            END DO
360         END DO
361      END DO
[825]362
[4161]363      ! Ice salinity, temperature and heat content
364      DO jk = 1, nlay_i
365         DO jl = 1, jpl ! loop over categories
366            DO jj = 1, jpj
367               DO ji = 1, jpi
[5123]368                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 
369                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
370                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
[825]371
[4161]372                   ! heat content per unit volume
[4688]373                   e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
[5123]374                      +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
375                      -   rcp     * ( ztmelts - rt0 ) )
[825]376
[5123]377                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
378                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
[5202]379               END DO
380            END DO
381         END DO
382      END DO
[825]383
[4688]384      tn_ice (:,:,:) = t_su (:,:,:)
385
386      ELSE 
[5123]387         ! if ln_iceini=false
[4688]388         a_i  (:,:,:) = 0._wp
389         v_i  (:,:,:) = 0._wp
390         v_s  (:,:,:) = 0._wp
391         smv_i(:,:,:) = 0._wp
392         oa_i (:,:,:) = 0._wp
393         ht_i (:,:,:) = 0._wp
394         ht_s (:,:,:) = 0._wp
395         sm_i (:,:,:) = 0._wp
396         o_i  (:,:,:) = 0._wp
397
398         e_i(:,:,:,:) = 0._wp
399         e_s(:,:,:,:) = 0._wp
400
401         DO jl = 1, jpl
402            DO jk = 1, nlay_i
[5123]403               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
[4688]404            END DO
405            DO jk = 1, nlay_s
[5123]406               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
[4688]407            END DO
408         END DO
409     
[5123]410      ENDIF ! ln_iceini
[4688]411     
412      at_i (:,:) = 0.0_wp
413      DO jl = 1, jpl
414         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
415      END DO
416      !
[825]417      !--------------------------------------------------------------------
[4161]418      ! 4) Global ice variables for output diagnostics                    |
[825]419      !--------------------------------------------------------------------
[4161]420      u_ice (:,:)     = 0._wp
421      v_ice (:,:)     = 0._wp
422      stress1_i(:,:)  = 0._wp
423      stress2_i(:,:)  = 0._wp
424      stress12_i(:,:) = 0._wp
[825]425
426      !--------------------------------------------------------------------
[4161]427      ! 5) Moments for advection
[825]428      !--------------------------------------------------------------------
429
[4161]430      sxopw (:,:) = 0._wp 
431      syopw (:,:) = 0._wp 
432      sxxopw(:,:) = 0._wp 
433      syyopw(:,:) = 0._wp 
434      sxyopw(:,:) = 0._wp
[3610]435
[4161]436      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
437      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
438      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
439      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
440      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
[825]441
[4161]442      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
443      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
444      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
445      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
446      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
[825]447
[4161]448      sxsal  (:,:,:)  = 0._wp
449      sysal  (:,:,:)  = 0._wp
450      sxxsal (:,:,:)  = 0._wp
451      syysal (:,:,:)  = 0._wp
452      sxysal (:,:,:)  = 0._wp
[825]453
[4335]454      sxage  (:,:,:)  = 0._wp
455      syage  (:,:,:)  = 0._wp
456      sxxage (:,:,:)  = 0._wp
457      syyage (:,:,:)  = 0._wp
458      sxyage (:,:,:)  = 0._wp
459
[825]460
[4688]461      CALL wrk_dealloc( jpi, jpj, zswitch )
[4161]462      CALL wrk_dealloc( jpi, jpj, zhemis )
[4688]463      CALL wrk_dealloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
464      CALL wrk_dealloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
[4161]465
[825]466   END SUBROUTINE lim_istate
467
468   SUBROUTINE lim_istate_init
469      !!-------------------------------------------------------------------
470      !!                   ***  ROUTINE lim_istate_init  ***
471      !!       
472      !! ** Purpose : Definition of initial state of the ice
473      !!
[4161]474      !! ** Method : Read the namiceini namelist and check the parameter
475      !!       values called at the first timestep (nit000)
[825]476      !!
[4161]477      !! ** input :
478      !!        Namelist namiceini
479      !!
480      !! history :
481      !!  8.5  ! 03-08 (C. Ethe) original code
482      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
[825]483      !!-----------------------------------------------------------------------------
[5123]484      NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s,  &
485         &                                      rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s
[4147]486      INTEGER :: ios                 ! Local integer output status for namelist read
[825]487      !!-----------------------------------------------------------------------------
[2528]488      !
[4147]489      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
490      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
491901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
492
493      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
494      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
495902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
[4624]496      IF(lwm) WRITE ( numoni, namiceini )
[4161]497
498      ! Define the initial parameters
499      ! -------------------------
500
501      IF(lwp) THEN
[825]502         WRITE(numout,*)
503         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
504         WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5123]505         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
506         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
507         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
508         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
509         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
510         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
511         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
512         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
513         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
514         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
515         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
516         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
[825]517      ENDIF
[4161]518
[825]519   END SUBROUTINE lim_istate_init
520
521#else
522   !!----------------------------------------------------------------------
523   !!   Default option :         Empty module          NO LIM sea-ice model
524   !!----------------------------------------------------------------------
525CONTAINS
526   SUBROUTINE lim_istate          ! Empty routine
527   END SUBROUTINE lim_istate
528#endif
529
530   !!======================================================================
531END MODULE limistate
Note: See TracBrowser for help on using the repository browser.