source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 6772

Last change on this file since 6772 was 6772, checked in by cbricaud, 4 years ago

clean in coarsening branch

  • Property svn:keywords set to Id
File size: 44.2 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   USE fldread          ! read input fields
31   USE iom
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) ::   rn_thres_sst   ! threshold water temperature for initial sea ice
40   REAL(wp) ::   rn_hts_ini_n   ! initial snow thickness in the north
41   REAL(wp) ::   rn_hts_ini_s   ! initial snow thickness in the south
42   REAL(wp) ::   rn_hti_ini_n   ! initial ice thickness in the north
43   REAL(wp) ::   rn_hti_ini_s   ! initial ice thickness in the south
44   REAL(wp) ::   rn_ati_ini_n   ! initial leads area in the north
45   REAL(wp) ::   rn_ati_ini_s   ! initial leads area in the south
46   REAL(wp) ::   rn_smi_ini_n   ! initial salinity
47   REAL(wp) ::   rn_smi_ini_s   ! initial salinity
48   REAL(wp) ::   rn_tmi_ini_n   ! initial temperature
49   REAL(wp) ::   rn_tmi_ini_s   ! initial temperature
50
51   INTEGER , PARAMETER ::   jpfldi    = 7           ! maximum number of files to read
52   INTEGER , PARAMETER ::   jp_hicif = 1           ! index of thick (m)    at T-point
53   INTEGER , PARAMETER ::   jp_hsnif = 2           ! index of thick (m)    at T-point
54   INTEGER , PARAMETER ::   jp_frld  = 3           ! index of ice fraction (%) at T-point
55   INTEGER , PARAMETER ::   jp_sist  = 4           ! index of ice surface temp (K)    at T-point
56   INTEGER , PARAMETER ::   jp_tbif1 = 5           ! index of ice temp lev1 (K) at T-point
57   INTEGER , PARAMETER ::   jp_tbif2 = 6           ! index of ice temp lev2 (K) at T-point
58   INTEGER , PARAMETER ::   jp_tbif3 = 7           ! index of ice temp lev3 (K) at T-point
59   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si    ! structure of input fields (file informations, fields read)
60
61   REAL(wp),DIMENSION(:,:)  ,ALLOCATABLE :: hicif_ini,hsnif_ini,frld_ini,sist_ini, zswitch
62   REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: tbif_ini
63
64   LOGICAL  ::  ln_iceini    ! initialization or not
65   LOGICAL  ::  ln_limini_file   ! Ice initialization state from 2D netcdf file
66   !!----------------------------------------------------------------------
67   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
68   !! $Id$
69   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   SUBROUTINE lim_istate
74      !!-------------------------------------------------------------------
75      !!                    ***  ROUTINE lim_istate  ***
76      !!
77      !! ** Purpose :   defined the sea-ice initial state
78      !!
79      !! ** Method  :   
80      !!                This routine will put some ice where ocean
81      !!                is at the freezing point, then fill in ice
82      !!                state variables using prescribed initial
83      !!                values in the namelist           
84      !!
85      !! ** Steps   :   
86      !!                1) Read namelist
87      !!                2) Basal temperature; ice and hemisphere masks
88      !!                3) Fill in the ice thickness distribution using gaussian
89      !!                4) Fill in space-dependent arrays for state variables
90      !!                5) Diagnostic arrays
91      !!                6) Lateral boundary conditions
92      !!
93      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
94      !!              where there is no ice (clem: I do not know why, is it mandatory?)
95      !!
96      !! History :
97      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
98      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
99      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
100      !!--------------------------------------------------------------------
101
102      !! * Local variables
103      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
104      REAL(wp)   :: ztmelts, zdh
105      INTEGER    :: i_hemis, i_fill, jl0 
106      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv
107      REAL(wp), POINTER, DIMENSION(:)     :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini
108      REAL(wp), POINTER, DIMENSION(:,:)   :: zh_i_ini, za_i_ini, zv_i_ini
109      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index
110      !--------------------------------------------------------------------
111
112      CALL wrk_alloc( jpi, jpj, zhemis )
113      CALL wrk_alloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
114      CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
115
116      IF(lwp) WRITE(numout,*)
117      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
118      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
119
120      !--------------------------------------------------------------------
121      ! 1) Read namelist
122      !--------------------------------------------------------------------
123
124      CALL lim_istate_init     !  reading the initials parameters of the ice
125
126      ! surface temperature
127      DO jl = 1, jpl ! loop over categories
128         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
129         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
130      END DO
131
132      ! basal temperature (considered at freezing point)
133      t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1) 
134
135      IF( ln_iceini ) THEN
136
137      !--------------------------------------------------------------------
138      ! 2) Basal temperature, ice mask and hemispheric index
139      !--------------------------------------------------------------------
140
141      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
142         DO ji = 1, jpi
143            IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN
144               zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
145            ELSE                                                                                   
146               zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
147            ENDIF
148         END DO
149      END DO
150
151
152      ! Hemispheric index
153      DO jj = 1, jpj
154         DO ji = 1, jpi
155            IF( fcor(ji,jj) >= 0._wp ) THEN   
156               zhemis(ji,jj) = 1 ! Northern hemisphere
157            ELSE
158               zhemis(ji,jj) = 2 ! Southern hemisphere
159            ENDIF
160         END DO
161      END DO
162
163      !--------------------------------------------------------------------
164      ! 3) Initialization of sea ice state variables
165      !--------------------------------------------------------------------
166      IF( ln_limini_file )THEN
167
168         CALL limini_file
169
170      ELSE
171
172      !-----------------------------
173      ! 3.1) Hemisphere-dependent arrays
174      !-----------------------------
175      ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
176      zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s  ! ice thickness
177      zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s  ! snow depth
178      zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s  ! ice concentration
179      zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s  ! bulk ice salinity
180      ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s  ! temperature (ice and snow)
181
182      zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:)   ! ice volume
183
184      !---------------------------------------------------------------------
185      ! 3.2) Distribute ice concentration and thickness into the categories
186      !---------------------------------------------------------------------
187      ! a gaussian distribution for ice concentration is used
188      ! then we check whether the distribution fullfills
189      ! volume and area conservation, positivity and ice categories bounds
190      DO i_hemis = 1, 2
191
192      ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
193
194      ! note for the great nemo engineers:
195      ! only very few of the WRITE statements are necessary for the reference version
196      ! they were one day useful, but now i personally doubt of their
197      ! potential for bringing anything useful
198
199      DO i_fill = jpl, 1, -1
200         IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
201            !----------------------------
202            ! fill the i_fill categories
203            !----------------------------
204            ! *** 1 category to fill
205            IF ( i_fill .EQ. 1 ) THEN
206               zh_i_ini(1,i_hemis)       = zht_i_ini(i_hemis)
207               za_i_ini(1,i_hemis)       = zat_i_ini(i_hemis)
208               zh_i_ini(2:jpl,i_hemis)   = 0._wp
209               za_i_ini(2:jpl,i_hemis)   = 0._wp
210            ELSE
211
212               ! *** >1 categores to fill
213               !--- Ice thicknesses in the i_fill - 1 first categories
214               DO jl = 1, i_fill - 1
215                  zh_i_ini(jl,i_hemis) = hi_mean(jl)
216               END DO
217               
218               !--- jl0: most likely index where cc will be maximum
219               DO jl = 1, jpl
220                  IF ( ( zht_i_ini(i_hemis) >  hi_max(jl-1) ) .AND. &
221                     & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN
222                     jl0 = jl
223                  ENDIF
224               END DO
225               jl0 = MIN(jl0, i_fill)
226               
227               !--- Concentrations
228               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl))
229               DO jl = 1, i_fill - 1
230                  IF ( jl .NE. jl0 ) THEN
231                     zsigma               = 0.5 * zht_i_ini(i_hemis)
232                     zarg                 = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma
233                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2)
234                  ENDIF
235               END DO
236               
237               zA = 0. ! sum of the areas in the jpl categories
238               DO jl = 1, i_fill - 1
239                 zA = zA + za_i_ini(jl,i_hemis)
240               END DO
241               za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category
242               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
243         
244               !--- Ice thickness in the last category
245               zV = 0. ! sum of the volumes of the N-1 categories
246               DO jl = 1, i_fill - 1
247                  zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis)
248               END DO
249               zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 
250               IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
251
252               !--- volumes
253               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis)
254               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
255
256            ENDIF ! i_fill
257
258            !---------------------
259            ! Compatibility tests
260            !---------------------
261            ! Test 1: area conservation
262            zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons )
263            IF ( zconv .LT. 1.0e-6 ) THEN
264               ztest_1 = 1
265            ELSE 
266              ! this write is useful
267              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 
268               ztest_1 = 0
269            ENDIF
270
271            ! Test 2: volume conservation
272            zV_cons = SUM(zv_i_ini(:,i_hemis))
273            zconv = ABS(zvt_i_ini(i_hemis) - zV_cons)
274
275            IF ( zconv .LT. 1.0e-6 ) THEN
276               ztest_2 = 1
277            ELSE
278              ! this write is useful
279              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
280                            ' zvt_i_ini = ', zvt_i_ini(i_hemis)
281               ztest_2 = 0
282            ENDIF
283
284            ! Test 3: thickness of the last category is in-bounds ?
285            IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN
286               ztest_3 = 1
287            ELSE
288               ! this write is useful
289               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', &
290               zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
291               ztest_3 = 0
292            ENDIF
293
294            ! Test 4: positivity of ice concentrations
295            ztest_4 = 1
296            DO jl = 1, jpl
297               IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 
298                  ! this write is useful
299                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis)
300                  ztest_4 = 0
301               ENDIF
302            END DO
303
304         ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
305 
306         ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
307
308      END DO ! i_fill
309
310      IF(lwp) THEN
311         WRITE(numout,*) ' ztests : ', ztests
312         IF ( ztests .NE. 4 ) THEN
313            WRITE(numout,*)
314            WRITE(numout,*) ' !!!! ALERT                  !!! '
315            WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
316            WRITE(numout,*)
317            WRITE(numout,*) ' *** ztests is not equal to 4 '
318            WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
319            WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis)
320            WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis)
321         ENDIF ! ztests .NE. 4
322      ENDIF
323     
324      END DO ! i_hemis
325
326      !---------------------------------------------------------------------
327      ! 3.3) Space-dependent arrays for ice state variables
328      !---------------------------------------------------------------------
329
330      ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
331      DO jl = 1, jpl ! loop over categories
332         DO jj = 1, jpj
333            DO ji = 1, jpi
334               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
335               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness
336               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth
337               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity
338               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day)
339               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
340
341               ! This case below should not be used if (ht_s/ht_i) is ok in namelist
342               ! In case snow load is in excess that would lead to transformation from snow to ice
343               ! Then, transfer the snow excess into the ice (different from limthd_dh)
344               zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
345               ! recompute ht_i, ht_s avoiding out of bounds values
346               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
347               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
348
349               ! ice volume, salt content, age content
350               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
351               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
352               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
353               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
354            END DO
355         END DO
356      END DO
357
358      ! Snow temperature and heat content
359      DO jk = 1, nlay_s
360         DO jl = 1, jpl ! loop over categories
361            DO jj = 1, jpj
362               DO ji = 1, jpi
363                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0
364                   ! Snow energy of melting
365                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
366
367                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
368                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
369               END DO
370            END DO
371         END DO
372      END DO
373
374      ! Ice salinity, temperature and heat content
375      DO jk = 1, nlay_i
376         DO jl = 1, jpl ! loop over categories
377            DO jj = 1, jpj
378               DO ji = 1, jpi
379                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 
380                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
381                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
382
383                   ! heat content per unit volume
384                   e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
385                      +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
386                      -   rcp     * ( ztmelts - rt0 ) )
387
388                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
389                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
390               END DO
391            END DO
392         END DO
393      END DO
394
395      tn_ice (:,:,:) = t_su (:,:,:)
396
397      ENDIF !ln_limini_file
398
399      ELSE 
400         ! if ln_iceini=false
401         a_i  (:,:,:) = 0._wp
402         v_i  (:,:,:) = 0._wp
403         v_s  (:,:,:) = 0._wp
404         smv_i(:,:,:) = 0._wp
405         oa_i (:,:,:) = 0._wp
406         ht_i (:,:,:) = 0._wp
407         ht_s (:,:,:) = 0._wp
408         sm_i (:,:,:) = 0._wp
409         o_i  (:,:,:) = 0._wp
410
411         e_i(:,:,:,:) = 0._wp
412         e_s(:,:,:,:) = 0._wp
413
414         DO jl = 1, jpl
415            DO jk = 1, nlay_i
416               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
417            END DO
418            DO jk = 1, nlay_s
419               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
420            END DO
421         END DO
422
423      ENDIF ! ln_iceini
424     
425      at_i (:,:) = 0.0_wp
426      DO jl = 1, jpl
427         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
428      END DO
429      !
430      !--------------------------------------------------------------------
431      ! 4) Global ice variables for output diagnostics                    |
432      !--------------------------------------------------------------------
433      u_ice (:,:)     = 0._wp
434      v_ice (:,:)     = 0._wp
435      stress1_i(:,:)  = 0._wp
436      stress2_i(:,:)  = 0._wp
437      stress12_i(:,:) = 0._wp
438
439      !--------------------------------------------------------------------
440      ! 5) Moments for advection
441      !--------------------------------------------------------------------
442
443      sxopw (:,:) = 0._wp 
444      syopw (:,:) = 0._wp 
445      sxxopw(:,:) = 0._wp 
446      syyopw(:,:) = 0._wp 
447      sxyopw(:,:) = 0._wp
448
449      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
450      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
451      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
452      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
453      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
454
455      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
456      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
457      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
458      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
459      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
460
461      sxsal  (:,:,:)  = 0._wp
462      sysal  (:,:,:)  = 0._wp
463      sxxsal (:,:,:)  = 0._wp
464      syysal (:,:,:)  = 0._wp
465      sxysal (:,:,:)  = 0._wp
466
467      sxage  (:,:,:)  = 0._wp
468      syage  (:,:,:)  = 0._wp
469      sxxage (:,:,:)  = 0._wp
470      syyage (:,:,:)  = 0._wp
471      sxyage (:,:,:)  = 0._wp
472
473
474      CALL wrk_dealloc( jpi, jpj, zhemis )
475      CALL wrk_dealloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini )
476      CALL wrk_dealloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
477
478   END SUBROUTINE lim_istate
479
480   SUBROUTINE lim_istate_init
481      !!-------------------------------------------------------------------
482      !!                   ***  ROUTINE lim_istate_init  ***
483      !!       
484      !! ** Purpose : Definition of initial state of the ice
485      !!
486      !! ** Method : Read the namiceini namelist and check the parameter
487      !!       values called at the first timestep (nit000)
488      !!
489      !! ** input :
490      !!        Namelist namiceini
491      !!
492      !! history :
493      !!  8.5  ! 03-08 (C. Ethe) original code
494      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
495      !!-----------------------------------------------------------------------------
496      !
497      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
498      INTEGER :: ji,jj
499      INTEGER :: ifpr, ierror
500      !
501      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ice files
502      TYPE(FLD_N)                    ::   sn_hicif, sn_hsnif, sn_frld, sn_sist
503      TYPE(FLD_N)                    ::   sn_tbif1, sn_tbif2, sn_tbif3
504      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
505      !
506      NAMELIST/namiceini/ ln_iceini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
507         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
508         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
509         &                sn_hicif, sn_hsnif, sn_frld, sn_sist,                                 &
510         &                sn_tbif1, sn_tbif2, sn_tbif3, cn_dir
511      !!-----------------------------------------------------------------------------
512      !
513      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
514      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
515901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
516
517      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
518      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
519902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
520      IF(lwm) WRITE ( numoni, namiceini )
521
522      slf_i(jp_hicif) = sn_hicif  ;  slf_i(jp_hsnif) = sn_hsnif
523      slf_i(jp_frld)  = sn_frld   ;  slf_i(jp_sist)  = sn_sist
524      slf_i(jp_tbif1) = sn_tbif1  ;  slf_i(jp_tbif2) = sn_tbif2  ; slf_i(jp_tbif3) = sn_tbif3
525
526      ! Define the initial parameters
527      ! -------------------------
528
529      IF(lwp) THEN
530         WRITE(numout,*)
531         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
532         WRITE(numout,*) '~~~~~~~~~~~~~~~'
533         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
534         WRITE(numout,*) '   initialization with ice (T) or not (F)   ln_limini_file  = ', ln_limini_file
535         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
536         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
537         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
538         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
539         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
540         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
541         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
542         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
543         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
544         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
545         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
546      ENDIF
547
548      IF( ln_limini_file ) THEN                      ! Ice initialization using input file
549         !
550         ierr = alloc_lim_istate_init()
551         !
552!         CALL iom_open( 'Ice_initialization.nc', inum_ice )
553!         !
554!         IF( inum_ice > 0 ) THEN
555!            IF(lwp) WRITE(numout,*)
556!            IF(lwp) WRITE(numout,*) '                  ice state initialization with : Ice_initialization.nc'
557!
558!            CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif_ini )
559!            CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif_ini )
560!            CALL iom_get( inum_ice, jpdom_data, 'frld' , frld_ini  )
561!            CALL iom_get( inum_ice, jpdom_data, 'ts'   , sist_ini  )
562!            CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif_ini(1:nlci,1:nlcj,:),   &
563!                 &        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,3 /) )
564!            ! put some values in the extra-halo...
565
566         ! set si structure
567         ALLOCATE( si(jpfldi), STAT=ierror )
568         IF( ierror > 0 ) THEN
569            CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' )   ;   RETURN
570         ENDIF
571
572         DO ifpr= 1, jpfldi
573            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
574            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
575         END DO
576
577         ! fill si with slf_i and control print
578         CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' )
579
580         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
581
582         hicif_ini(:,:)  = si(jp_hicif)%fnow(:,:,1)
583         hsnif_ini(:,:)  = si(jp_hsnif)%fnow(:,:,1)
584         frld_ini(:,:)   = si(jp_frld)%fnow(:,:,1)
585         sist_ini(:,:)   = si(jp_sist)%fnow(:,:,1)
586         tbif_ini(:,:,1) = si(jp_tbif1)%fnow(:,:,1)
587         tbif_ini(:,:,2) = si(jp_tbif2)%fnow(:,:,1)
588         tbif_ini(:,:,3) = si(jp_tbif3)%fnow(:,:,1)
589
590         DO jj = nlcj+1, jpj   ;   tbif_ini(1:nlci,jj,:) = tbif_ini(1:nlci,nlej,:)   ;   END DO
591         DO ji = nlci+1, jpi   ;   tbif_ini(ji    ,: ,:) = tbif_ini(nlei  ,:   ,:)   ;   END DO
592
593!            CALL iom_close( inum_ice)
594!            !
595!         ENDIF
596      ENDIF
597
598   END SUBROUTINE lim_istate_init
599
600   SUBROUTINE limini_file
601      !!-----------------------------------------------------------------------------
602      !!
603      !!
604      !!
605      !!
606      !!-----------------------------------------------------------------------------
607      INTEGER    :: jl,ji,jj,jk
608      INTEGER    :: jl0
609      INTEGER    :: i_fill,jit,jjt
610      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv,zH
611      REAL(wp)   :: eps=1.e-6
612      REAL(wp)   :: zmin,zmax
613      !rbb REAL(wp)   :: epsi20,ztmelts,zdh
614      REAL(wp)   ::ztmelts,zdh
615
616      REAL(wp), POINTER, DIMENSION(:,:)   :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini
617      REAL(wp), POINTER, DIMENSION(:,:,:) :: zv_i_ini
618      REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ini,za_i_ini
619      REAL(wp), POINTER, DIMENSION(:,:)   :: zidto    ! ice indicator
620       !-----------------------------------------------------------------------------
621      IF(lwp)WRITE(numout,*)"limistate: read file : "
622
623      CALL wrk_alloc(jpl,jpi,jpj, zv_i_ini)
624      CALL wrk_alloc(    jpi,jpj, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini )
625      CALL wrk_alloc(    jpl,jpi,jpj,zht_i_ini,za_i_ini)
626      CALL wrk_alloc(    jpi,jpj,zidto )
627
628      zhm_i_ini(:,:) = hicif_ini(:,:)  ! ice thickness
629      zat_i_ini(:,:) = 1._wp - frld_ini(:,:)   ! ice concentration
630      zvt_i_ini(:,:) = zhm_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume
631      zhm_s_ini(:,:) = hsnif_ini(:,:)  ! snow depth
632
633      zht_i_ini(:,:,:) = 0._wp
634      za_i_ini(:,:,:) = 0._wp
635      zv_i_ini(:,:,:) = 0._wp
636
637      zat_i_ini(:,:) = MIN( zat_i_ini(:,:) , 1.0_wp )
638
639
640      DO ji = 1, jpi
641      DO jj = 1, jpj
642
643         IF( zat_i_ini(ji,jj) .GT. 0._wp .AND. zhm_i_ini(ji,jj) .GT. 0._wp )THEN
644
645
646            IF( gphit(ji,jj) .GE. 0._wp )THEN ; zsm_i_ini(ji,jj) = rn_smi_ini_n
647            ELSE                              ; zsm_i_ini(ji,jj) = rn_smi_ini_s
648            ENDIF
649
650            jl0 = 1
651            DO jl = 2, jpl
652               IF ( ( zhm_i_ini(ji,jj) .GT. hi_max(jl-1) ) .AND. &
653                  (   zhm_i_ini(ji,jj) .LE. hi_max(jl)   )       ) THEN
654               jl0 = jl
655               ENDIF
656            END DO
657
658            IF( jl0==1 )THEN
659
660               zht_i_ini(1,ji,jj)       = zhm_i_ini(ji,jj)
661               za_i_ini(1,ji,jj)        = zat_i_ini(ji,jj)
662               zht_i_ini(2:jpl,ji,jj)   = 0._wp
663               za_i_ini(2:jpl,ji,jj)    = 0._wp
664
665            ELSE ! jl0 ne 1
666               ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
667
668               DO i_fill = jpl, 1, -1
669                  IF( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
670
671                     !----------------------------
672                     ! fill the i_fill categories
673                     !----------------------------
674                     ! *** 1 category to fill
675                     IF( i_fill .EQ. 1 ) THEN
676                        zht_i_ini(1,ji,jj)       = zhm_i_ini(ji,jj)
677                        za_i_ini(1,ji,jj)        = zat_i_ini(ji,jj)
678                        zht_i_ini(2:jpl,ji,jj)   = 0._wp
679                        za_i_ini(2:jpl,ji,jj)    = 0._wp
680                     ELSE
681
682                        ! *** >1 categores to fill
683                        !--- Ice thicknesses in the i_fill - 1 first categories
684                        DO jl = 1, i_fill - 1
685                           zht_i_ini(jl,ji,jj)    = 0.5 * ( hi_max(jl) + hi_max(jl-1) )
686                        END DO
687
688                        !--- jl0: most likely index where cc will be maximum
689                        DO jl = 1, jpl
690                        IF ( ( zhm_i_ini(ji,jj) .GT. hi_max(jl-1) ) .AND. &
691                              ( zhm_i_ini(ji,jj) .LE. hi_max(jl)   ) ) THEN
692                            jl0 = jl
693                        ENDIF
694                        END DO
695                        jl0 = MIN(jl0, i_fill)
696
697                        !--- Concentrations
698                        za_i_ini(jl0,ji,jj)      = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
699                        DO jl = 1, i_fill - 1
700                        IF ( jl .NE. jl0 ) THEN
701                             zsigma               = 0.5 * zhm_i_ini(ji,jj)
702                             zarg                 = ( zht_i_ini(jl,ji,jj) - zhm_i_ini(ji,jj) ) / zsigma
703                             za_i_ini(jl,ji,jj) = za_i_ini(jl0,ji,jj) * EXP(-zarg**2)
704                        ENDIF
705                        END DO
706
707                        zA = 0. ! sum of the areas in the jpl categories
708                        DO jl = 1, i_fill - 1
709                           zA = zA + za_i_ini(jl,ji,jj)
710                        END DO
711                        za_i_ini(i_fill,ji,jj)   = zat_i_ini(ji,jj) - zA ! ice conc in the last category
712                        IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, ji,jj) = 0._wp
713
714                        !--- Ice thickness in the last category
715                        zV = 0. ! sum of the volumes of the N-1 categories
716                        DO jl = 1, i_fill - 1
717                           zV = zV + za_i_ini(jl,ji,jj)*zht_i_ini(jl,ji,jj)
718                        END DO
719                        zht_i_ini(i_fill,ji,jj) = ( zvt_i_ini(ji,jj) - zV ) /za_i_ini(i_fill,ji,jj)
720                        IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, ji,jj) = 0._wp
721
722                        !--- volumes
723                        zv_i_ini(:,ji,jj) = za_i_ini(:,ji,jj) * zht_i_ini(:,ji,jj)
724                        IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, ji,jj) = 0._wp
725
726                     ENDIF ! i_fill
727
728                     !---------------------
729                     ! Compatibility tests
730                     !---------------------
731                     ! Test 1: area conservation
732                     zA_cons = SUM(za_i_ini(:,ji,jj)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons )
733                     IF ( zconv .LT. 1.0e-6 ) THEN
734                        ztest_1 = 1
735                     ELSE
736                      ! this write is useful
737                      !WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj)
738                      !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
739                      !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
740                      !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
741                      !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
742                      !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
743                      !WRITE(numout,*) ' hi_max ',hi_max
744                      !WRITE(numout,*) ' jl0 = ',jl0
745                      !WRITE(numout,*) ' vol = ',zvt_i_ini(ji,jj),SUM(zv_i_ini(:,ji,jj))
746                      ztest_1 = 0
747                     ENDIF
748
749                     ! Test 2: volume conservation
750                     zV_cons = SUM(zv_i_ini(:,ji,jj))
751                     zconv = ABS(zvt_i_ini(ji,jj) - zV_cons)
752
753                     IF ( zconv .LT. 1.0e-6 ) THEN
754                        ztest_2 = 1
755                     ELSE
756                        ! this write is useful
757                        !WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
758                        !    ' zvt_i_ini = ', zvt_i_ini(ji,jj)
759                        !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
760                        !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
761                        !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
762                        !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
763                        !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
764                        !WRITE(numout,*) ' hi_max ',hi_max
765                        !WRITE(numout,*) ' jl0 = ',jl0
766                        ztest_2 = 0
767                     ENDIF
768
769                     ! Test 3: thickness of the last category is in-bounds ?
770                     IF ( zht_i_ini(i_fill, ji,jj) .GT. hi_max(i_fill-1) ) THEN
771                     ztest_3 = 1
772                     ELSE
773                     ! this write is useful
774                     !WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,ji,jj) = ', &
775                     !zht_i_ini(i_fill,ji,jj), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
776                     !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
777                     !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
778                     !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
779                     !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
780                     !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
781                     !WRITE(numout,*) ' hi_max ',hi_max
782                     !WRITE(numout,*) ' jl0 = ',jl0
783                     ztest_3 = 0
784                     ENDIF
785
786                     ! Test 4: positivity of ice concentrations
787                     ztest_4 = 1
788                     DO jl = 1, jpl
789                     IF ( za_i_ini(jl,ji,jj) .LT. 0._wp ) THEN
790                        ! this write is useful
791                        !WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, 'WITH A = ', za_i_ini(jl,ji,jj)
792                        !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
793                        !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
794                        !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
795                        !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
796                        !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
797                        !WRITE(numout,*) ' hi_max ',hi_max
798                        !WRITE(numout,*) ' jl0 = ',jl0
799                        !WRITE(numout,*)
800                        ztest_4 = 0
801                     ENDIF
802                     END DO
803
804                  ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
805
806                  ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
807
808               END DO ! i_fill
809
810               !WRITE(numout,*) ' ztests : ', ztests
811               !IF ( ztests .NE. 4 ) THEN
812               !WRITE(numout,*)
813               !WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
814               !WRITE(numout,*) ' !!!! RED ALERT                  !!! '
815               !WRITE(numout,*) ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!'
816               !WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
817               !WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
818               !WRITE(numout,*) 'ji,jj,narea ',ji,jj,narea
819               !WRITE(numout,*) ' *** ztests is not equal to 4 '
820               !WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2,ztest_3,ztest_4
821               !WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
822               !WRITE(numout,*) ' zhm_i_ini : ', zhm_i_ini(ji,jj)
823               !WRITE(numout,*) ' zht_i_ini(:,jij,jj) ',zht_i_ini(:,ji,jj)
824               !WRITE(numout,*) ' za_i_ini(:,jij,jj) ',za_i_ini(:,ji,jj)
825               !WRITE(numout,*) ' hi_max ',hi_max
826               !ENDIF ! ztests .NE. 4
827
828            ENDIF  !  jl0 ne 1
829
830         ENDIF  !  zat_i_ini ne 0
831      END DO ! jj
832      END DO ! ji
833
834
835      !---------------------------------------------------------------------
836      ! 3.3) Space-dependent arrays for ice state variables
837      !---------------------------------------------------------------------
838
839      ! Ice concentration, thickness and volume, ice salinity, ice age, surface
840      ! temperature
841      DO jl = 1, jpl ! loop over categories
842         DO jj = 1, jpj
843            DO ji = 1, jpi
844               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,ji,jj)  ! concentration
845               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zht_i_ini(jl,ji,jj)   !ice thickness
846
847               IF( zhm_i_ini( ji,jj ) .GT. 0_wp )THEN ; ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zhm_s_ini( ji,jj ) / zhm_i_ini( ji,jj ) )
848               ELSE                                   ; ht_s(ji,jj,jl)  = 0._wp
849               ENDIF
850               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj) !+ (1._wp - zswitch(ji,jj) ) * rn_simin ! salinity
851               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp + ( 1._wp -zswitch(ji,jj) ) ! age
852               t_su(ji,jj,jl)  = sist_ini(ji,jj)
853
854               ! This case below should not be used if (ht_s/ht_i) is ok in
855               ! namelist
856               ! In case snow load is in excess that would lead to
857               ! transformation from snow to ice
858               ! Then, transfer the snow excess into the ice (different from
859               ! limthd_dh)
860               zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) *ht_i(ji,jj,jl) ) * r1_rau0 )
861               ! recompute ht_i, ht_s avoiding out of bounds values
862               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
863               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic *r1_rhosn )
864
865               ! ice volume, salt content, age content
866               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              !ice volume
867               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              !snow volume
868               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) *v_i(ji,jj,jl) ! salt content
869               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               !age content
870            END DO ! ji
871         END DO ! jj
872      END DO ! jl
873
874      !cbr
875      DO jk = 1, nlay_s
876         DO  jl = 1, jpl ! loop over categories
877            !rbb t_s(:,:,1,jl) =  tbif_ini(:,:,1)
878            t_s(:,:,1,jl) =  tbif_ini(:,:,1)*zswitch(:,:)+ ( 1._wp - zswitch(:,:) ) * rt0
879         END DO ! jl
880      END DO ! jk
881
882      ! Snow temperature and heat content
883      DO jk = 1, nlay_s
884         DO jl = 1, jpl ! loop over categories
885            DO jj = 1, jpj
886               DO ji = 1, jpi
887!cbr???                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
888                   ! Snow energy of melting
889                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
890
891                   ! Mutliply by volume, and divide by number of layers to get
892                   ! heat content in J/m2
893                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) *r1_nlay_s
894               END DO ! ji
895            END DO ! jj
896         END DO ! jl
897      END DO ! jk
898
899      ! Ice salinity, temperature and heat content
900      DO  jk = 1, nlay_i
901         DO jl = 1, jpl ! loop over categories
902            DO jj = 1, jpj
903               DO ji = 1, jpi
904!cbr???                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
905                   t_i(ji,jj,jk,jl) =  tbif_ini(ji,jj,2)*zswitch(ji,jj)+ ( 1._wp - zswitch(ji,jj) ) * rt0
906                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
907                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0           !Melting temperature in K
908
909                   ! heat content per unit volume
910                   e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
911                      +   lfus    * ( 1._wp - (ztmelts-rt0) /MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
912                      -   rcp     * ( ztmelts - rt0 ) )
913
914                   ! Mutliply by ice volume, and divide by number of layers to
915                   ! get heat content in J/m2
916                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
917               END DO ! ji
918            END DO ! jj
919         END DO ! jl
920      END DO ! jk
921
922      !cbr tmp CALL wrk_dealloc(jpl,jpi,jpj, zht_i_ini, za_i_ini, zv_i_ini)
923      CALL wrk_dealloc(jpl,jpi,jpj, zv_i_ini)
924      CALL wrk_dealloc(    jpl,jpi,jpj,zht_i_ini,za_i_ini)
925      CALL wrk_dealloc(    jpi,jpj, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini,zsm_i_ini )
926      CALL wrk_dealloc(    jpi,jpj,zidto )
927
928  END SUBROUTINE limini_file
929
930
931  INTEGER FUNCTION alloc_lim_istate_init()
932      !!-----------------------------------------------------------------------------
933      !!
934      !!
935      !!
936      !!
937      !!-----------------------------------------------------------------------------
938      INTEGER :: ierr(1)
939      !!-----------------------------------------------------------------------------
940      ALLOCATE( hicif_ini(jpi,jpj) , hsnif_ini(jpi,jpj) , frld_ini(jpi,jpj) , sist_ini(jpi,jpj) , zswitch(jpi,jpj) , tbif_ini(jpi,jpj,3) , Stat=ierr(1) )
941      alloc_lim_istate_init = MAXVAL(ierr)
942      IF( alloc_lim_istate_init /= 0 )   CALL ctl_warn( 'lim_istate_init: failed to allocate arrays')
943
944   END FUNCTION alloc_lim_istate_init
945#else
946   !!----------------------------------------------------------------------
947   !!   Default option :         Empty module          NO LIM sea-ice model
948   !!----------------------------------------------------------------------
949CONTAINS
950   SUBROUTINE lim_istate          ! Empty routine
951   END SUBROUTINE lim_istate
952#endif
953
954   !!======================================================================
955END MODULE limistate
Note: See TracBrowser for help on using the repository browser.