source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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