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 NEMO/branches/UKMO/r8395_restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: NEMO/branches/UKMO/r8395_restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 10758

Last change on this file since 10758 was 10758, checked in by jcastill, 5 years ago

Remove svn keywords

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