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.
iceistate.F90 in branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90 @ 8752

Last change on this file since 8752 was 8752, checked in by dancopsey, 7 years ago

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

File size: 28.7 KB
Line 
1MODULE iceistate
2   !!======================================================================
3   !!                     ***  MODULE  iceistate  ***
4   !!   sea-ice : Initialization of 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'                                       ESIM 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, ln_ice_embd 
23   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
24   USE eosbn2         ! equation of state
25   USE domvvl         ! Variable volume
26   USE ice            ! sea-ice variables
27   USE icevar         ! ice_var_salprof
28   !
29   USE in_out_manager ! I/O manager
30   USE iom            ! I/O manager library
31   USE lib_mpp        ! MPP library
32   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
33   USE fldread        ! read input fields
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   ice_istate        ! called by icestp.F90
39   PUBLIC   ice_istate_init   ! called by icestp.F90
40
41   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
42   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point
43   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point
44   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
45   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
46   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
47   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
48   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
49   !
50   ! ** namelist (namini) **
51   LOGICAL  ::   ln_iceini        ! initialization or not
52   LOGICAL  ::   ln_iceini_file   ! Ice initialization state from 2D netcdf file
53   REAL(wp) ::   rn_thres_sst     ! threshold water temperature for initial sea ice
54   REAL(wp) ::   rn_hts_ini_n     ! initial snow thickness in the north
55   REAL(wp) ::   rn_hts_ini_s     ! initial snow thickness in the south
56   REAL(wp) ::   rn_hti_ini_n     ! initial ice thickness in the north
57   REAL(wp) ::   rn_hti_ini_s     ! initial ice thickness in the south
58   REAL(wp) ::   rn_ati_ini_n     ! initial leads area in the north
59   REAL(wp) ::   rn_ati_ini_s     ! initial leads area in the south
60   REAL(wp) ::   rn_smi_ini_n     ! initial salinity
61   REAL(wp) ::   rn_smi_ini_s     ! initial salinity
62   REAL(wp) ::   rn_tmi_ini_n     ! initial temperature
63   REAL(wp) ::   rn_tmi_ini_s     ! initial temperature
64   
65   !!----------------------------------------------------------------------
66   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
67   !! $Id: iceistate.F90 8378 2017-07-26 13:55:59Z clem $
68   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
69   !!----------------------------------------------------------------------
70CONTAINS
71
72   SUBROUTINE ice_istate
73      !!-------------------------------------------------------------------
74      !!                    ***  ROUTINE ice_istate  ***
75      !!
76      !! ** Purpose :   defined the sea-ice initial state
77      !!
78      !! ** Method  :   This routine will put some ice where ocean
79      !!                is at the freezing point, then fill in ice
80      !!                state variables using prescribed initial
81      !!                values in the namelist           
82      !!
83      !! ** Steps   :   1) Set initial surface and basal temperatures
84      !!                2) Recompute or read sea ice state variables
85      !!                3) Fill in the ice thickness distribution using gaussian
86      !!                4) Fill in space-dependent arrays for state variables
87      !!                5) snow-ice mass computation
88      !!                6) store before fields
89      !!
90      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
91      !!              where there is no ice (clem: I do not know why, is it mandatory?)
92      !!--------------------------------------------------------------------
93      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
94      REAL(wp) ::   ztmelts, zdh
95      INTEGER  ::   i_hemis, i_fill, jl0
96      REAL(wp) ::   zarg, zV, zconv, zdv, zfac
97      INTEGER , DIMENSION(4)           ::   itest
98      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
99      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
100      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
101      REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
102      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini, za_i_ini                         !data by cattegories to fill
103      !--------------------------------------------------------------------
104
105      IF(lwp) WRITE(numout,*)
106      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
107      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
108
109      !--------------------------------------------------------------------
110      ! 1) Set surface and bottom temperatures to initial values
111      !--------------------------------------------------------------------
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      IF( ln_iceini ) THEN
124         !-----------------------------------------------------------
125         ! 2) Compute or read sea ice variables ===> single category
126         !-----------------------------------------------------------
127         !
128         !                             !---------------!
129         IF( ln_iceini_file )THEN      ! Read a file   !
130            !                          !---------------!
131            !
132            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
133            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
134            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
135            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
136            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
137            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
138            !
139            WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 
140            ELSEWHERE                       ; zswitch(:,:) = 0._wp
141            END WHERE
142            zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)
143            !
144         !                             !---------------!
145         ELSE                          ! Read namelist !
146            !                          !---------------!
147
148           ! no ice if sst <= t-freez + ttest
149            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 
150            ELSEWHERE                                                                  ; zswitch(:,:) = tmask(:,:,1)
151            END WHERE
152
153            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
154            WHERE( ff_t(:,:) >= 0._wp )
155               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
156               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
157               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
158               zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
159               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
160               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
161            ELSEWHERE
162               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
163               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
164               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
165               zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
166               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
167               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
168            END WHERE
169            zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)
170            !
171         ENDIF
172         
173         !------------------------------------------------------------------
174         ! 3) Distribute ice concentration and thickness into the categories
175         !------------------------------------------------------------------
176         ! a gaussian distribution for ice concentration is used
177         ! then we check whether the distribution fullfills
178         ! volume and area conservation, positivity and ice categories bounds
179         zh_i_ini(:,:,:) = 0._wp 
180         za_i_ini(:,:,:) = 0._wp
181         !
182         DO jj = 1, jpj
183            DO ji = 1, jpi
184               !
185               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
186
187                  ! find which category (jl0) the input ice thickness falls into
188                  jl0 = jpl
189                  DO jl = 1, jpl
190                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
191                        jl0 = jl
192                        CYCLE
193                     ENDIF
194                  END DO
195                  !
196                  itest(:) = 0
197                  i_fill   = jpl + 1                                            !------------------------------------
198                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
199                     !                                                          !------------------------------------
200                     i_fill = i_fill - 1
201                     !
202                     zh_i_ini(ji,jj,:) = 0._wp 
203                     za_i_ini(ji,jj,:) = 0._wp
204                     itest(:) = 0
205                     !
206                     IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
207                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
208                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
209                     ELSE                         !-- case ice is thicker: fill categories >1
210                        ! thickness
211                        DO jl = 1, i_fill-1
212                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
213                        END DO
214                        !
215                        ! concentration
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                        ! last 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                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
227                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
228
229                        ! clem: correction if concentration of upper cat is greater than lower cat
230                        !       (it should be a gaussian around jl0 but sometimes it is not)
231                        IF ( jl0 /= jpl ) THEN
232                           DO jl = jpl, jl0+1, -1
233                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
234                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
235                                 zh_i_ini(ji,jj,jl    ) = 0._wp
236                                 za_i_ini(ji,jj,jl    ) = 0._wp
237                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
238                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
239                              END IF
240                           ENDDO
241                        ENDIF
242                        !
243                     ENDIF
244
245                     ! Compatibility tests
246                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation
247                     IF ( zconv < epsi06 ) itest(1) = 1
248                     
249                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation
250                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
251                     IF ( zconv < epsi06 ) itest(2) = 1
252                     
253                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ?
254                     
255                     itest(4) = 1
256                     DO jl = 1, i_fill
257                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations
258                     END DO
259                     !                                                          !----------------------------
260                  END DO                                                        ! end iteration on categories
261                  !                                                             !----------------------------
262                  !
263                  IF( lwp .AND. SUM(itest) /= 4 ) THEN
264                     WRITE(numout,*)
265                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
266                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
267                     WRITE(numout,*)
268                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
269                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
270                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
271                  ENDIF
272               
273               ENDIF
274               !
275            END DO   
276         END DO   
277
278         !---------------------------------------------------------------------
279         ! 4) Fill in sea ice arrays
280         !---------------------------------------------------------------------
281
282         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
283         DO jl = 1, jpl ! loop over categories
284            DO jj = 1, jpj
285               DO ji = 1, jpi
286                  a_i(ji,jj,jl)  = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
287                  h_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
288                  s_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
289                  o_i(ji,jj,jl)  = 0._wp                                                     ! age (0 day)
290                  t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
291
292                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
293                    h_s(ji,jj,jl)= h_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
294                  ELSE
295                    h_s(ji,jj,jl)= 0._wp
296                  ENDIF
297
298                  ! This case below should not be used if (h_s/h_i) is ok in namelist
299                  ! In case snow load is in excess that would lead to transformation from snow to ice
300                  ! Then, transfer the snow excess into the ice (different from icethd_dh)
301                  zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 
302                  ! recompute h_i, h_s avoiding out of bounds values
303                  h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
304                  h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
305
306                  ! ice volume, salt content, age content
307                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
308                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
309                  sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
310                  oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
311               END DO
312            END DO
313         END DO
314
315         ! for constant salinity in time
316         IF( nn_icesal /= 2 )  THEN
317            CALL ice_var_salprof
318            sv_i = s_i * v_i
319         ENDIF
320           
321         ! Snow temperature and heat content
322         DO jk = 1, nlay_s
323            DO jl = 1, jpl ! loop over categories
324               DO jj = 1, jpj
325                  DO ji = 1, jpi
326                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
327                     ! Snow energy of melting
328                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
329
330                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
331                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
332                  END DO
333               END DO
334            END DO
335         END DO
336
337         ! Ice salinity, temperature and heat content
338         DO jk = 1, nlay_i
339            DO jl = 1, jpl ! loop over categories
340               DO jj = 1, jpj
341                  DO ji = 1, jpi
342                     t_i (ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
343                     sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
344                     ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
345
346                     ! heat content per unit volume
347                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
348                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
349                        -   rcp     * ( ztmelts - rt0 ) )
350
351                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
352                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
353                  END DO
354               END DO
355            END DO
356         END DO
357
358         tn_ice (:,:,:) = t_su (:,:,:)
359
360         ! Melt pond volume and fraction
361         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp
362         ELSE                                     ;   zfac = 0._wp   ;   ENDIF
363         DO jl = 1, jpl
364            a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac
365            h_ip     (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac
366         END DO
367         a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:) 
368         v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:)
369
370      ELSE ! if ln_iceini=false
371         a_i  (:,:,:) = 0._wp
372         v_i  (:,:,:) = 0._wp
373         v_s  (:,:,:) = 0._wp
374         sv_i (:,:,:) = 0._wp
375         oa_i (:,:,:) = 0._wp
376         h_i  (:,:,:) = 0._wp
377         h_s  (:,:,:) = 0._wp
378         s_i  (:,:,:) = 0._wp
379         o_i  (:,:,:) = 0._wp
380
381         e_i(:,:,:,:) = 0._wp
382         e_s(:,:,:,:) = 0._wp
383
384         DO jl = 1, jpl
385            DO jk = 1, nlay_i
386               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
387            END DO
388            DO jk = 1, nlay_s
389               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
390            END DO
391         END DO
392
393         a_ip(:,:,:)      = 0._wp
394         v_ip(:,:,:)      = 0._wp
395         a_ip_frac(:,:,:) = 0._wp
396         h_ip     (:,:,:) = 0._wp
397
398      ENDIF ! ln_iceini
399     
400      at_i (:,:) = 0.0_wp
401      DO jl = 1, jpl
402         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
403      END DO
404      !
405      ! --- set ice velocities --- !
406      u_ice (:,:)     = 0._wp
407      v_ice (:,:)     = 0._wp
408      !
409      !----------------------------------------------
410      ! 5) Snow-ice mass (case ice is fully embedded)
411      !----------------------------------------------
412      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhosn * v_s(:,:,:) + rhoic * v_i(:,:,:), dim=3  )   ! snow+ice mass
413      snwice_mass_b(:,:) = snwice_mass(:,:)
414      !
415      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
416
417         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
418         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
419
420         IF( .NOT.ln_linssh ) THEN
421           
422            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
423            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
424         
425            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
426               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
427               e3t_b(:,:,jk) = e3t_n(:,:,jk)
428               e3t_a(:,:,jk) = e3t_n(:,:,jk)
429            END DO
430           
431            ! Reconstruction of all vertical scale factors at now and before time-steps
432            ! =========================================================================
433            ! Horizontal scale factor interpolations
434            ! --------------------------------------
435            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
436            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
437            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
438            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
439            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
440            ! Vertical scale factor interpolations
441            ! ------------------------------------
442            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
443            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
444            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
445            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
446            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
447            ! t- and w- points depth
448            ! ----------------------
449            !!gm not sure of that....
450            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
451            gdepw_n(:,:,1) = 0.0_wp
452            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
453            DO jk = 2, jpk
454               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
455               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
456               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
457            END DO
458         ENDIF
459      ENDIF
460     
461      !------------------------------------
462      ! 6) store fields at before time-step
463      !------------------------------------
464      ! it is only necessary for the 1st interpolation by Agrif
465      a_i_b  (:,:,:)   = a_i  (:,:,:)
466      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
467      v_i_b  (:,:,:)   = v_i  (:,:,:)
468      v_s_b  (:,:,:)   = v_s  (:,:,:)
469      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
470      sv_i_b (:,:,:)   = sv_i (:,:,:)
471      oa_i_b (:,:,:)   = oa_i (:,:,:)
472      u_ice_b(:,:)     = u_ice(:,:)
473      v_ice_b(:,:)     = v_ice(:,:)
474
475!!!clem
476!!      ! Output the initial state and forcings
477!!      CALL dia_wri_state( 'output.init', nit000 )
478!!!     
479
480   END SUBROUTINE ice_istate
481
482   SUBROUTINE ice_istate_init
483      !!-------------------------------------------------------------------
484      !!                   ***  ROUTINE ice_istate_init  ***
485      !!       
486      !! ** Purpose : Definition of initial state of the ice
487      !!
488      !! ** Method : Read the namini namelist and check the parameter
489      !!       values called at the first timestep (nit000)
490      !!
491      !! ** input :
492      !!        Namelist namini
493      !!
494      !! history :
495      !!  8.5  ! 03-08 (C. Ethe) original code
496      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
497      !!-----------------------------------------------------------------------------
498      !
499      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
500      INTEGER :: ji,jj
501      INTEGER :: ifpr, ierror
502      !
503      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
504      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
505      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
506      !
507      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
508         &             rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
509         &             rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
510         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
511      !!-----------------------------------------------------------------------------
512      !
513      REWIND( numnam_ice_ref )              ! Namelist namini in reference namelist : Ice initial state
514      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
515901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist', lwp )
516
517      REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state
518      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
519902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist', lwp )
520      IF(lwm) WRITE ( numoni, namini )
521
522      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
523      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
524      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
525      !
526      !
527      IF(lwp) THEN                          ! control print
528         WRITE(numout,*)
529         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
530         WRITE(numout,*) '~~~~~~~~~~~~~~~'
531         WRITE(numout,*) '   Namelist namini:'
532         WRITE(numout,*) '      initialization with ice (T) or not (F)                 ln_iceini       = ', ln_iceini
533         WRITE(numout,*) '      ice initialization from a netcdf file                  ln_iceini_file  = ', ln_iceini_file
534         WRITE(numout,*) '      max delta ocean temp. above Tfreeze with initial ice   rn_thres_sst    = ', rn_thres_sst
535         WRITE(numout,*) '      initial snow thickness in the north                    rn_hts_ini_n    = ', rn_hts_ini_n
536         WRITE(numout,*) '      initial snow thickness in the south                    rn_hts_ini_s    = ', rn_hts_ini_s 
537         WRITE(numout,*) '      initial ice thickness  in the north                    rn_hti_ini_n    = ', rn_hti_ini_n
538         WRITE(numout,*) '      initial ice thickness  in the south                    rn_hti_ini_s    = ', rn_hti_ini_s
539         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_n    = ', rn_ati_ini_n
540         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_s    = ', rn_ati_ini_s
541         WRITE(numout,*) '      initial  ice salinity  in the north                    rn_smi_ini_n    = ', rn_smi_ini_n
542         WRITE(numout,*) '      initial  ice salinity  in the south                    rn_smi_ini_s    = ', rn_smi_ini_s
543         WRITE(numout,*) '      initial  ice/snw temp  in the north                    rn_tmi_ini_n    = ', rn_tmi_ini_n
544         WRITE(numout,*) '      initial  ice/snw temp  in the south                    rn_tmi_ini_s    = ', rn_tmi_ini_s
545      ENDIF
546
547      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
548         !
549         ! set si structure
550         ALLOCATE( si(jpfldi), STAT=ierror )
551         IF( ierror > 0 ) THEN
552            CALL ctl_stop( 'Ice_ini in iceistate: unable to allocate si structure' )   ;   RETURN
553         ENDIF
554
555         DO ifpr = 1, jpfldi
556            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
557            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
558         END DO
559
560         ! fill si with slf_i and control print
561         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' )
562
563         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
564
565      ENDIF
566
567   END SUBROUTINE ice_istate_init
568
569#else
570   !!----------------------------------------------------------------------
571   !!   Default option :         Empty module         NO ESIM sea-ice model
572   !!----------------------------------------------------------------------
573#endif
574
575   !!======================================================================
576END MODULE iceistate
Note: See TracBrowser for help on using the repository browser.