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

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

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