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/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 4900

Last change on this file since 4900 was 4900, checked in by cetlod, 10 years ago

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