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

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

implement several developments for LIM3: new advection scheme (ultimate-macho, not yet perfect) ; lateral ice melt ; enabling/disabling thermo and dyn with namelist options ; simplifications (including a clarified namelist)
  • Property svn:keywords set to Id
File size: 28.0 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                          !this write is useful
264                          IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj) 
265                          ztest_1 = 0
266                        ENDIF
267
268                        ! Test 2: volume conservation
269                        zV_cons = SUM(zv_i_ini(ji,jj,:))
270                        zconv = ABS(zvt_i_ini(ji,jj) - zV_cons)
271
272                        IF( zconv .LT. 1.0e-6 ) THEN
273                           ztest_2 = 1
274                        ELSE
275                           !this write is useful
276                           IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
277                                                    ' zvt_i_ini = ', zvt_i_ini(ji,jj)
278                           ztest_2 = 0
279                        ENDIF
280
281                        ! Test 3: thickness of the last category is in-bounds ?
282                        IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN
283                           ztest_3 = 1
284                        ELSE
285                           ! this write is useful
286                           IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(ji,jj,i_fill) = ', &
287                           zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
288                           IF(lwp) WRITE(numout,*) ' ji,jj,i_fill ',ji,jj,i_fill
289                           IF(lwp) WRITE(numout,*) 'zht_i_ini ',zht_i_ini(ji,jj)
290                           ztest_3 = 0
291                        ENDIF
292
293                        ! Test 4: positivity of ice concentrations
294                        ztest_4 = 1
295                        DO jl = 1, jpl
296                           IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN 
297                              ! this write is useful
298                              IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(ji,jj,jl)
299                              ztest_4 = 0
300                           ENDIF
301                        END DO
302
303                     ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
304 
305                     ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
306
307                  END DO ! i_fill
308
309                  IF(lwp) THEN
310                     WRITE(numout,*) ' ztests : ', ztests
311                     IF( ztests .NE. 4 )THEN
312                        WRITE(numout,*)
313                        WRITE(numout,*) ' !!!! ALERT                  !!! '
314                        WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
315                        WRITE(numout,*)
316                        WRITE(numout,*) ' *** ztests is not equal to 4 '
317                        WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
318                        WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
319                        WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
320                     ENDIF ! ztests .NE. 4
321                  ENDIF
322     
323               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp
324
325            ENDDO   
326         ENDDO   
327
328         !---------------------------------------------------------------------
329         ! 3.3) Space-dependent arrays for ice state variables
330         !---------------------------------------------------------------------
331
332         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
333         DO jl = 1, jpl ! loop over categories
334            DO jj = 1, jpj
335               DO ji = 1, jpi
336                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
337                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
338                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
339                  o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                                    ! age (1 day)
340                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
341
342                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
343                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
344                  ELSE
345                    ht_s(ji,jj,jl)= 0._wp
346                  ENDIF
347
348                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist
349                  ! In case snow load is in excess that would lead to transformation from snow to ice
350                  ! Then, transfer the snow excess into the ice (different from limthd_dh)
351                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
352                  ! recompute ht_i, ht_s avoiding out of bounds values
353                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
354                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
355
356                  ! ice volume, salt content, age content
357                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
358                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
359                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
360                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
361               END DO
362            END DO
363         END DO
364
365         ! for constant salinity in time
366         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
367            CALL lim_var_salprof
368            smv_i = sm_i * v_i
369         ENDIF
370           
371         ! Snow temperature and heat content
372         DO jk = 1, nlay_s
373            DO jl = 1, jpl ! loop over categories
374               DO jj = 1, jpj
375                  DO ji = 1, jpi
376                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
377                     ! Snow energy of melting
378                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
379
380                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
381                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
382                  END DO
383               END DO
384            END DO
385         END DO
386
387         ! Ice salinity, temperature and heat content
388         DO jk = 1, nlay_i
389            DO jl = 1, jpl ! loop over categories
390               DO jj = 1, jpj
391                  DO ji = 1, jpi
392                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
393                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
394                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
395
396                     ! heat content per unit volume
397                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
398                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
399                        -   rcp     * ( ztmelts - rt0 ) )
400
401                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
402                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
403                  END DO
404               END DO
405            END DO
406         END DO
407
408         tn_ice (:,:,:) = t_su (:,:,:)
409
410      ELSE ! if ln_limini=false
411         a_i  (:,:,:) = 0._wp
412         v_i  (:,:,:) = 0._wp
413         v_s  (:,:,:) = 0._wp
414         smv_i(:,:,:) = 0._wp
415         oa_i (:,:,:) = 0._wp
416         ht_i (:,:,:) = 0._wp
417         ht_s (:,:,:) = 0._wp
418         sm_i (:,:,:) = 0._wp
419         o_i  (:,:,:) = 0._wp
420
421         e_i(:,:,:,:) = 0._wp
422         e_s(:,:,:,:) = 0._wp
423
424         DO jl = 1, jpl
425            DO jk = 1, nlay_i
426               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
427            END DO
428            DO jk = 1, nlay_s
429               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
430            END DO
431         END DO
432
433      ENDIF ! ln_limini
434     
435      at_i (:,:) = 0.0_wp
436      DO jl = 1, jpl
437         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
438      END DO
439      !
440      !--------------------------------------------------------------------
441      ! 4) Global ice variables for output diagnostics                    |
442      !--------------------------------------------------------------------
443      u_ice (:,:)     = 0._wp
444      v_ice (:,:)     = 0._wp
445      stress1_i(:,:)  = 0._wp
446      stress2_i(:,:)  = 0._wp
447      stress12_i(:,:) = 0._wp
448
449      !--------------------------------------------------------------------
450      ! 5) Moments for advection
451      !--------------------------------------------------------------------
452
453      sxopw (:,:) = 0._wp 
454      syopw (:,:) = 0._wp 
455      sxxopw(:,:) = 0._wp 
456      syyopw(:,:) = 0._wp 
457      sxyopw(:,:) = 0._wp
458
459      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
460      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
461      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
462      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
463      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
464
465      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
466      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
467      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
468      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
469      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
470
471      sxsal  (:,:,:)  = 0._wp
472      sysal  (:,:,:)  = 0._wp
473      sxxsal (:,:,:)  = 0._wp
474      syysal (:,:,:)  = 0._wp
475      sxysal (:,:,:)  = 0._wp
476
477      sxage  (:,:,:)  = 0._wp
478      syage  (:,:,:)  = 0._wp
479      sxxage (:,:,:)  = 0._wp
480      syyage (:,:,:)  = 0._wp
481      sxyage (:,:,:)  = 0._wp
482
483
484      CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini )
485      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 )
486      CALL wrk_dealloc( jpi, jpj,      zswitch )
487
488   END SUBROUTINE lim_istate
489
490   SUBROUTINE lim_istate_init
491      !!-------------------------------------------------------------------
492      !!                   ***  ROUTINE lim_istate_init  ***
493      !!       
494      !! ** Purpose : Definition of initial state of the ice
495      !!
496      !! ** Method : Read the namiceini namelist and check the parameter
497      !!       values called at the first timestep (nit000)
498      !!
499      !! ** input :
500      !!        Namelist namiceini
501      !!
502      !! history :
503      !!  8.5  ! 03-08 (C. Ethe) original code
504      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
505      !!-----------------------------------------------------------------------------
506      !
507      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
508      INTEGER :: ji,jj
509      INTEGER :: ifpr, ierror
510      !
511      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ice files
512      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
513      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
514      !
515      NAMELIST/namiceini/ ln_limini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
516         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
517         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
518         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
519      !!-----------------------------------------------------------------------------
520      !
521      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
522      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
523901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
524
525      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
526      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
527902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
528      IF(lwm) WRITE ( numoni, namiceini )
529
530      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
531      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
532      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
533
534      ! Define the initial parameters
535      ! -------------------------
536
537      IF(lwp) THEN
538         WRITE(numout,*)
539         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
540         WRITE(numout,*) '~~~~~~~~~~~~~~~'
541         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_limini     = ', ln_limini
542         WRITE(numout,*) '   ice initialization from a netcdf file      ln_limini_file  = ', ln_limini_file
543         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
544         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
545         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
546         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
547         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
548         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
549         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
550         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
551         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
552         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
553         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
554      ENDIF
555
556      IF( ln_limini_file ) THEN                      ! Ice initialization using input file
557         !
558         ! set si structure
559         ALLOCATE( si(jpfldi), STAT=ierror )
560         IF( ierror > 0 ) THEN
561            CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' )   ;   RETURN
562         ENDIF
563
564         DO ifpr = 1, jpfldi
565            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
566            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
567         END DO
568
569         ! fill si with slf_i and control print
570         CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' )
571
572         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
573
574      ENDIF
575
576   END SUBROUTINE lim_istate_init
577
578#else
579   !!----------------------------------------------------------------------
580   !!   Default option :         Empty module          NO LIM sea-ice model
581   !!----------------------------------------------------------------------
582CONTAINS
583   SUBROUTINE lim_istate          ! Empty routine
584   END SUBROUTINE lim_istate
585#endif
586
587   !!======================================================================
588END MODULE limistate
Note: See TracBrowser for help on using the repository browser.