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

source: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 5208

Last change on this file since 5208 was 5208, checked in by davestorkey, 9 years ago

Merge in changes from trunk up to 5021.

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