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

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 4166

Last change on this file since 4166 was 4166, checked in by cetlod, 11 years ago

dev_LOCEAN_2013 : bugs corrections, see ticket #1169

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