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

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 4924

Last change on this file since 4924 was 4924, checked in by mathiot, 9 years ago

UKM02_ice_shelves merged and SETTE tested with revision 4879 of trunk

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