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 @ 6746

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

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