source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 5541

Last change on this file since 5541 was 5541, checked in by jchanut, 5 years ago

change eos_fzp function into a subroutine (needed for AGRIF), #1530

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