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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 4528

Last change on this file since 4528 was 4335, checked in by clem, 10 years ago

debug: intialization of ice age moments for advection

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