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

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 6853

Last change on this file since 6853 was 6853, checked in by clem, 8 years ago

correct bugs on LIM3 rheology and outputs. Make sure the ice thickness distribution (used in initialization and bdy boundary conditions) is a gaussian.

  • Property svn:keywords set to Id
File size: 26.9 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 limvar           ! lim_var_salprof
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!!!clem
34!!   USE diawri
35!!!
36   
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   lim_istate      ! routine called by lim_init.F90
41
42   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
43   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point
44   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point
45   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
46   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
47   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
48   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
49   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
50   !!----------------------------------------------------------------------
51   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
52   !! $Id$
53   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE lim_istate
58      !!-------------------------------------------------------------------
59      !!                    ***  ROUTINE lim_istate  ***
60      !!
61      !! ** Purpose :   defined the sea-ice initial state
62      !!
63      !! ** Method  :   
64      !!                This routine will put some ice where ocean
65      !!                is at the freezing point, then fill in ice
66      !!                state variables using prescribed initial
67      !!                values in the namelist           
68      !!
69      !! ** Steps   :   
70      !!                1) Read namelist
71      !!                2) Basal temperature; ice and hemisphere masks
72      !!                3) Fill in the ice thickness distribution using gaussian
73      !!                4) Fill in space-dependent arrays for state variables
74      !!                5) Diagnostic arrays
75      !!                6) Lateral boundary conditions
76      !!
77      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
78      !!              where there is no ice (clem: I do not know why, is it mandatory?)
79      !!
80      !! History :
81      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
82      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
83      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
84      !!--------------------------------------------------------------------
85
86      !! * Local variables
87      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
88      REAL(wp)   :: ztmelts, zdh
89      INTEGER    :: i_hemis, i_fill, jl0 
90      REAL(wp)   :: zarg, zV, zconv, zdv 
91      REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator
92      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
93      REAL(wp), POINTER, DIMENSION(:,:)   :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
94      REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini                         !data by cattegories to fill
95      INTEGER , POINTER, DIMENSION(:)     :: itest
96      !--------------------------------------------------------------------
97
98      CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini )
99      CALL wrk_alloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
100      CALL wrk_alloc( jpi, jpj,      zswitch )
101      Call wrk_alloc( 4,             itest )
102
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      ! surface temperature
114      DO jl = 1, jpl ! loop over categories
115         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
116         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
117      END DO
118
119      ! basal temperature (considered at freezing point)
120      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
121      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
122
123
124      IF( ln_limini ) THEN
125
126         !--------------------------------------------------------------------
127         ! 2) Basal temperature, ice mask and hemispheric index
128         !--------------------------------------------------------------------
129
130         DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
131            DO ji = 1, jpi
132               IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst * tmask(ji,jj,1) ) THEN
133                  zswitch(ji,jj) = 0._wp                     ! no ice
134               ELSE                                                                                   
135                  zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
136               ENDIF
137            END DO
138         END DO
139
140         !--------------------------------------------------------------------
141         ! 3) Initialization of sea ice state variables
142         !--------------------------------------------------------------------
143         IF( ln_limini_file )THEN
144
145            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
146            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
147            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
148            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
149            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
150            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
151
152         ELSE ! ln_limini_file = F
153
154            !-----------------------------
155            ! 3.1) Hemisphere-dependent arrays
156            !-----------------------------
157            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
158            DO jj = 1, jpj
159               DO ji = 1, jpi
160                  IF( ff(ji,jj) >= 0._wp ) THEN
161                     zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj)
162                     zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj)
163                     zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj)
164                     zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
165                     zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj)
166                     ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
167                  ELSE
168                     zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj)
169                     zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj)
170                     zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj)
171                     zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
172                     zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj)
173                     ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
174                  ENDIF
175               END DO
176            END DO
177
178         ENDIF ! ln_limini_file
179         
180         zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume
181         !---------------------------------------------------------------------
182         ! 3.2) Distribute ice concentration and thickness into the categories
183         !---------------------------------------------------------------------
184         ! a gaussian distribution for ice concentration is used
185         ! then we check whether the distribution fullfills
186         ! volume and area conservation, positivity and ice categories bounds
187         zh_i_ini(:,:,:) = 0._wp 
188         za_i_ini(:,:,:) = 0._wp
189
190         DO jj = 1, jpj
191            DO ji = 1, jpi
192
193               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
194
195                  !--- jl0: most likely index where cc will be maximum
196                  jl0 = jpl
197                  DO jl = 1, jpl
198                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
199                        jl0 = jl
200                        CYCLE
201                     ENDIF
202                  END DO
203
204                  ! initialisation of tests
205                  itest(:)  = 0
206                 
207                  i_fill = jpl + 1                                             !====================================
208                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
209                     ! iteration                                               !====================================
210                     i_fill = i_fill - 1
211
212                     ! initialisation of ice variables for each try
213                     zh_i_ini(ji,jj,:) = 0._wp 
214                     za_i_ini(ji,jj,:) = 0._wp
215                     itest(:) = 0
216
217                     ! *** case very thin ice: fill only category 1
218                     IF ( i_fill == 1 ) THEN
219                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
220                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
221
222                     ! *** case ice is thicker: fill categories >1
223                     ELSE
224
225                        ! Fill ice thicknesses in the (i_fill-1) cat by hmean
226                        DO jl = 1, i_fill-1
227                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
228                        END DO
229               
230                        !--- Concentrations
231                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
232                        DO jl = 1, i_fill - 1
233                           IF( jl /= jl0 )THEN
234                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
235                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
236                           ENDIF
237                        END DO
238               
239                        ! Concentration in the last (i_fill) category
240                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
241
242                        ! Ice thickness in the last (i_fill) category
243                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
244                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
245
246                        ! clem: correction if concentration of upper cat is greater than lower cat
247                        !       (it should be a gaussian around jl0 but sometimes it is not)
248                        IF ( jl0 /= jpl ) THEN
249                           DO jl = jpl, jl0+1, -1
250                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
251                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
252                                 zh_i_ini(ji,jj,jl    ) = 0._wp
253                                 za_i_ini(ji,jj,jl    ) = 0._wp
254                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
255                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
256                              END IF
257                           ENDDO
258                        ENDIF
259
260                     ENDIF ! case ice is thick or thin
261
262                     !---------------------
263                     ! Compatibility tests
264                     !---------------------
265                     ! Test 1: area conservation
266                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )
267                     IF ( zconv < epsi06 ) itest(1) = 1
268                     
269                     ! Test 2: volume conservation
270                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &
271                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
272                     IF ( zconv < epsi06 ) itest(2) = 1
273                     
274                     ! Test 3: thickness of the last category is in-bounds ?
275                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
276                     
277                     ! Test 4: positivity of ice concentrations
278                     itest(4) = 1
279                     DO jl = 1, i_fill
280                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0
281                     END DO
282                     !                                      !============================
283                  END DO                                    ! end iteration on categories
284                  !                                         !============================
285
286                  IF( lwp .AND. SUM(itest) /= 4 ) THEN
287                     WRITE(numout,*)
288                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
289                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
290                     WRITE(numout,*)
291                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
292                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
293                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
294                  ENDIF
295               
296               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp
297           
298            ENDDO
299         ENDDO
300
301         !---------------------------------------------------------------------
302         ! 3.3) Space-dependent arrays for ice state variables
303         !---------------------------------------------------------------------
304
305         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
306         DO jl = 1, jpl ! loop over categories
307            DO jj = 1, jpj
308               DO ji = 1, jpi
309                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
310                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
311                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
312                  o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                                    ! age (1 day)
313                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
314
315                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
316                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
317                  ELSE
318                    ht_s(ji,jj,jl)= 0._wp
319                  ENDIF
320
321                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist
322                  ! In case snow load is in excess that would lead to transformation from snow to ice
323                  ! Then, transfer the snow excess into the ice (different from limthd_dh)
324                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
325                  ! recompute ht_i, ht_s avoiding out of bounds values
326                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
327                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
328
329                  ! ice volume, salt content, age content
330                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
331                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
332                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
333                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
334               END DO
335            END DO
336         END DO
337
338         ! for constant salinity in time
339         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
340            CALL lim_var_salprof
341            smv_i = sm_i * v_i
342         ENDIF
343           
344         ! Snow temperature and heat content
345         DO jk = 1, nlay_s
346            DO jl = 1, jpl ! loop over categories
347               DO jj = 1, jpj
348                  DO ji = 1, jpi
349                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
350                     ! Snow energy of melting
351                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
352
353                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
354                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
355                  END DO
356               END DO
357            END DO
358         END DO
359
360         ! Ice salinity, temperature and heat content
361         DO jk = 1, nlay_i
362            DO jl = 1, jpl ! loop over categories
363               DO jj = 1, jpj
364                  DO ji = 1, jpi
365                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
366                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
367                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
368
369                     ! heat content per unit volume
370                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
371                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
372                        -   rcp     * ( ztmelts - rt0 ) )
373
374                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
375                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
376                  END DO
377               END DO
378            END DO
379         END DO
380
381         tn_ice (:,:,:) = t_su (:,:,:)
382
383      ELSE ! if ln_limini=false
384         a_i  (:,:,:) = 0._wp
385         v_i  (:,:,:) = 0._wp
386         v_s  (:,:,:) = 0._wp
387         smv_i(:,:,:) = 0._wp
388         oa_i (:,:,:) = 0._wp
389         ht_i (:,:,:) = 0._wp
390         ht_s (:,:,:) = 0._wp
391         sm_i (:,:,:) = 0._wp
392         o_i  (:,:,:) = 0._wp
393
394         e_i(:,:,:,:) = 0._wp
395         e_s(:,:,:,:) = 0._wp
396
397         DO jl = 1, jpl
398            DO jk = 1, nlay_i
399               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
400            END DO
401            DO jk = 1, nlay_s
402               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
403            END DO
404         END DO
405
406      ENDIF ! ln_limini
407     
408      at_i (:,:) = 0.0_wp
409      DO jl = 1, jpl
410         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
411      END DO
412      !
413      !--------------------------------------------------------------------
414      ! 4) Global ice variables for output diagnostics                    |
415      !--------------------------------------------------------------------
416      u_ice (:,:)     = 0._wp
417      v_ice (:,:)     = 0._wp
418      stress1_i(:,:)  = 0._wp
419      stress2_i(:,:)  = 0._wp
420      stress12_i(:,:) = 0._wp
421
422      !--------------------------------------------------------------------
423      ! 5) Moments for advection
424      !--------------------------------------------------------------------
425
426      sxopw (:,:) = 0._wp 
427      syopw (:,:) = 0._wp 
428      sxxopw(:,:) = 0._wp 
429      syyopw(:,:) = 0._wp 
430      sxyopw(:,:) = 0._wp
431
432      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
433      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
434      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
435      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
436      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
437
438      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
439      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
440      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
441      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
442      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
443
444      sxsal  (:,:,:)  = 0._wp
445      sysal  (:,:,:)  = 0._wp
446      sxxsal (:,:,:)  = 0._wp
447      syysal (:,:,:)  = 0._wp
448      sxysal (:,:,:)  = 0._wp
449
450      sxage  (:,:,:)  = 0._wp
451      syage  (:,:,:)  = 0._wp
452      sxxage (:,:,:)  = 0._wp
453      syyage (:,:,:)  = 0._wp
454      sxyage (:,:,:)  = 0._wp
455
456!!!clem
457!!      ! Output the initial state and forcings
458!!      CALL dia_wri_state( 'output.init', nit000 )
459!!!
460     
461      CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini )
462      CALL wrk_dealloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
463      CALL wrk_dealloc( jpi, jpj,      zswitch )
464      Call wrk_dealloc( 4,             itest )
465
466   END SUBROUTINE lim_istate
467
468   SUBROUTINE lim_istate_init
469      !!-------------------------------------------------------------------
470      !!                   ***  ROUTINE lim_istate_init  ***
471      !!       
472      !! ** Purpose : Definition of initial state of the ice
473      !!
474      !! ** Method : Read the namiceini namelist and check the parameter
475      !!       values called at the first timestep (nit000)
476      !!
477      !! ** input :
478      !!        Namelist namiceini
479      !!
480      !! history :
481      !!  8.5  ! 03-08 (C. Ethe) original code
482      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
483      !!-----------------------------------------------------------------------------
484      !
485      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
486      INTEGER :: ji,jj
487      INTEGER :: ifpr, ierror
488      !
489      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ice files
490      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
491      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
492      !
493      NAMELIST/namiceini/ ln_limini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
494         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
495         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
496         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
497      !!-----------------------------------------------------------------------------
498      !
499      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
500      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
501901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
502
503      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
504      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
505902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
506      IF(lwm) WRITE ( numoni, namiceini )
507
508      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
509      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
510      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
511
512      ! Define the initial parameters
513      ! -------------------------
514
515      IF(lwp) THEN
516         WRITE(numout,*)
517         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
518         WRITE(numout,*) '~~~~~~~~~~~~~~~'
519         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_limini     = ', ln_limini
520         WRITE(numout,*) '   ice initialization from a netcdf file      ln_limini_file  = ', ln_limini_file
521         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
522         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
523         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
524         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
525         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
526         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
527         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
528         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
529         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
530         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
531         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
532      ENDIF
533
534      IF( ln_limini_file ) THEN                      ! Ice initialization using input file
535         !
536         ! set si structure
537         ALLOCATE( si(jpfldi), STAT=ierror )
538         IF( ierror > 0 ) THEN
539            CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' )   ;   RETURN
540         ENDIF
541
542         DO ifpr = 1, jpfldi
543            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
544            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
545         END DO
546
547         ! fill si with slf_i and control print
548         CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' )
549
550         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
551
552      ENDIF
553
554   END SUBROUTINE lim_istate_init
555
556#else
557   !!----------------------------------------------------------------------
558   !!   Default option :         Empty module          NO LIM sea-ice model
559   !!----------------------------------------------------------------------
560CONTAINS
561   SUBROUTINE lim_istate          ! Empty routine
562   END SUBROUTINE lim_istate
563#endif
564
565   !!======================================================================
566END MODULE limistate
Note: See TracBrowser for help on using the repository browser.