source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90 @ 8514

Last change on this file since 8514 was 8514, checked in by clem, 3 years ago

changes in style - part5 - almost done

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