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/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90 @ 9119

Last change on this file since 9119 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch

File size: 29.5 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      INTEGER  ::   i_hemis, i_fill, jl0   ! local integers
95      REAL(wp) ::   ztmelts, zdh
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)   [Kelvin]
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            ! 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            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
153            WHERE( ff_t(:,:) >= 0._wp )
154               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
155               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
156               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
157               zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
158               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
159               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
160            ELSEWHERE
161               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
162               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
163               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
164               zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
165               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
166               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
167            END WHERE
168            zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)
169            !
170         ENDIF
171         
172         !------------------------------------------------------------------
173         ! 3) Distribute ice concentration and thickness into the categories
174         !------------------------------------------------------------------
175         ! a gaussian distribution for ice concentration is used
176         ! then we check whether the distribution fullfills
177         ! volume and area conservation, positivity and ice categories bounds
178         zh_i_ini(:,:,:) = 0._wp 
179         za_i_ini(:,:,:) = 0._wp
180         !
181         DO jj = 1, jpj
182            DO ji = 1, jpi
183               !
184               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
185
186                  ! find which category (jl0) the input ice thickness falls into
187                  jl0 = jpl
188                  DO jl = 1, jpl
189                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
190                        jl0 = jl
191                        CYCLE
192                     ENDIF
193                  END DO
194                  !
195                  itest(:) = 0
196                  i_fill   = jpl + 1                                            !------------------------------------
197                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
198                     !                                                          !------------------------------------
199                     i_fill = i_fill - 1
200                     !
201                     zh_i_ini(ji,jj,:) = 0._wp 
202                     za_i_ini(ji,jj,:) = 0._wp
203                     itest(:) = 0
204                     !
205                     IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
206                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
207                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
208                     ELSE                         !-- case ice is thicker: fill categories >1
209                        ! thickness
210                        DO jl = 1, i_fill-1
211                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
212                        END DO
213                        !
214                        ! concentration
215                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
216                        DO jl = 1, i_fill - 1
217                           IF( jl /= jl0 )THEN
218                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
219                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
220                           ENDIF
221                        END DO
222
223                        ! last category
224                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
225                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
226                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
227
228                        ! correction if concentration of upper cat is greater than lower cat
229                        !   (it should be a gaussian around jl0 but sometimes it is not)
230                        IF ( jl0 /= jpl ) THEN
231                           DO jl = jpl, jl0+1, -1
232                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
233                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
234                                 zh_i_ini(ji,jj,jl    ) = 0._wp
235                                 za_i_ini(ji,jj,jl    ) = 0._wp
236                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
237                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
238                              END IF
239                           ENDDO
240                        ENDIF
241                        !
242                     ENDIF
243                     !
244                     ! Compatibility tests
245                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation
246                     IF ( zconv < epsi06 ) itest(1) = 1
247                     !
248                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation
249                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
250                     IF ( zconv < epsi06 ) itest(2) = 1
251                     !
252                     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 ?
253                     !
254                     itest(4) = 1
255                     DO jl = 1, i_fill
256                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations
257                     END DO
258                     !                                                          !----------------------------
259                  END DO                                                        ! end iteration on categories
260                  !                                                             !----------------------------
261                  IF( lwp .AND. SUM(itest) /= 4 ) THEN
262                     WRITE(numout,*)
263                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
264                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
265                     WRITE(numout,*)
266                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
267                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
268                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
269                  ENDIF
270                  !
271               ENDIF
272               !
273            END DO   
274         END DO   
275
276         !---------------------------------------------------------------------
277         ! 4) Fill in sea ice arrays
278         !---------------------------------------------------------------------
279         !
280         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
281         DO jl = 1, jpl ! loop over categories
282            DO jj = 1, jpj
283               DO ji = 1, jpi
284                  a_i(ji,jj,jl)  = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
285                  h_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
286                  s_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
287                  o_i(ji,jj,jl)  = 0._wp                                                     ! age (0 day)
288                  t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
289                  !
290                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
291                    h_s(ji,jj,jl)= h_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
292                  ELSE
293                    h_s(ji,jj,jl)= 0._wp
294                  ENDIF
295                  !
296                  ! This case below should not be used if (h_s/h_i) is ok in namelist
297                  ! In case snow load is in excess that would lead to transformation from snow to ice
298                  ! Then, transfer the snow excess into the ice (different from icethd_dh)
299                  zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 
300                  ! recompute h_i, h_s avoiding out of bounds values
301                  h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
302                  h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
303                  !
304                  ! ice volume, salt content, age content
305                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
306                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
307                  sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
308                  oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
309               END DO
310            END DO
311         END DO
312         !
313         IF( nn_icesal /= 2 )  THEN         ! for constant salinity in time
314            CALL ice_var_salprof
315            sv_i = s_i * v_i
316         ENDIF
317         
318         ! Snow temperature and heat content
319         DO jk = 1, nlay_s
320            DO jl = 1, jpl ! loop over categories
321               DO jj = 1, jpj
322                  DO ji = 1, jpi
323                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
324                     ! Snow energy of melting
325                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
326                     !
327                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
328                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
329                  END DO
330               END DO
331            END DO
332         END DO
333         !
334         ! Ice salinity, temperature and heat content
335         DO jk = 1, nlay_i
336            DO jl = 1, jpl ! loop over categories
337               DO jj = 1, jpj
338                  DO ji = 1, jpi
339                     t_i (ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
340                     sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
341                     ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
342                     !
343                     ! heat content per unit volume
344                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) )         &
345                        &             + lfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   &
346                        &             - rcp  * ( ztmelts - rt0 ) )
347                     !
348                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
349                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
350                  END DO
351               END DO
352            END DO
353         END DO
354         !
355         tn_ice (:,:,:) = t_su (:,:,:)
356         t1_ice (:,:,:) = t_i (:,:,1,:)   ! initialisation of 1st layer temp for coupled simu
357         cnd_ice(:,:,:) = 0._wp           ! initialisation of the effective conductivity at the top of ice/snow (Jules coupling)
358
359         ! Melt pond volume and fraction
360         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp
361         ELSE                                     ;   zfac = 0._wp
362         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         tn_ice (:,:,:) = t_i (:,:,1,:)
394         t1_ice (:,:,:) = t_i (:,:,1,:)   ! initialisation of 1st layer temp for coupled simu
395         cnd_ice(:,:,:) = 0._wp           ! initialisation of the effective conductivity at the top of ice/snow (Jules coupling)
396         
397         a_ip(:,:,:)      = 0._wp
398         v_ip(:,:,:)      = 0._wp
399         a_ip_frac(:,:,:) = 0._wp
400         h_ip     (:,:,:) = 0._wp
401         !
402      ENDIF ! ln_iceini
403      !
404      at_i (:,:) = 0.0_wp
405      DO jl = 1, jpl
406         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
407      END DO
408      !
409      ! --- set ice velocities --- !
410      u_ice (:,:) = 0._wp
411      v_ice (:,:) = 0._wp
412      !
413      !----------------------------------------------
414      ! 5) Snow-ice mass (case ice is fully embedded)
415      !----------------------------------------------
416      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhosn * v_s(:,:,:) + rhoic * v_i(:,:,:), dim=3  )   ! snow+ice mass
417      snwice_mass_b(:,:) = snwice_mass(:,:)
418      !
419      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
420         !
421         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
422         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
423         !
424         IF( .NOT.ln_linssh ) THEN
425            !
426            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
427            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
428            !
429            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
430               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
431               e3t_b(:,:,jk) = e3t_n(:,:,jk)
432               e3t_a(:,:,jk) = e3t_n(:,:,jk)
433            END DO
434            !
435            ! Reconstruction of all vertical scale factors at now and before time-steps
436            ! =========================================================================
437            ! Horizontal scale factor interpolations
438            ! --------------------------------------
439            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
440            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
441            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
442            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
443            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
444            ! Vertical scale factor interpolations
445            ! ------------------------------------
446            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
447            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
448            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
449            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
450            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
451            ! t- and w- points depth
452            ! ----------------------
453            !!gm not sure of that....
454            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
455            gdepw_n(:,:,1) = 0.0_wp
456            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
457            DO jk = 2, jpk
458               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
459               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
460               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
461            END DO
462         ENDIF
463      ENDIF
464     
465      !------------------------------------
466      ! 6) store fields at before time-step
467      !------------------------------------
468      ! it is only necessary for the 1st interpolation by Agrif
469      a_i_b  (:,:,:)   = a_i  (:,:,:)
470      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
471      v_i_b  (:,:,:)   = v_i  (:,:,:)
472      v_s_b  (:,:,:)   = v_s  (:,:,:)
473      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
474      sv_i_b (:,:,:)   = sv_i (:,:,:)
475      oa_i_b (:,:,:)   = oa_i (:,:,:)
476      u_ice_b(:,:)     = u_ice(:,:)
477      v_ice_b(:,:)     = v_ice(:,:)
478
479!!clem: output of initial state should be written here but it is impossible because
480!!      the ocean and ice are in the same file
481!!      CALL dia_wri_state( 'output.init', nit000 )
482      !
483   END SUBROUTINE ice_istate
484
485
486   SUBROUTINE ice_istate_init
487      !!-------------------------------------------------------------------
488      !!                   ***  ROUTINE ice_istate_init  ***
489      !!       
490      !! ** Purpose :   Definition of initial state of the ice
491      !!
492      !! ** Method  :   Read the namini namelist and check the parameter
493      !!              values called at the first timestep (nit000)
494      !!
495      !! ** input   :  Namelist namini
496      !!
497      !!-----------------------------------------------------------------------------
498      INTEGER ::   ji, jj
499      INTEGER ::   ios, ierr, inum_ice   ! Local integer output status for namelist read
500      INTEGER ::   ifpr, ierror
501      !
502      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
503      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
504      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
505      !
506      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
507         &             rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
508         &             rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
509         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
510      !!-----------------------------------------------------------------------------
511      !
512      REWIND( numnam_ice_ref )              ! Namelist namini in reference namelist : Ice initial state
513      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
514901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist', lwp )
515      !
516      REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state
517      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
518902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist', lwp )
519      IF(lwm) WRITE ( numoni, namini )
520      !
521      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
522      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
523      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
524      !
525      !
526      IF(lwp) THEN                          ! control print
527         WRITE(numout,*)
528         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
529         WRITE(numout,*) '~~~~~~~~~~~~~~~'
530         WRITE(numout,*) '   Namelist namini:'
531         WRITE(numout,*) '      initialization with ice (T) or not (F)                 ln_iceini       = ', ln_iceini
532         WRITE(numout,*) '      ice initialization from a netcdf file                  ln_iceini_file  = ', ln_iceini_file
533         WRITE(numout,*) '      max delta ocean temp. above Tfreeze with initial ice   rn_thres_sst    = ', rn_thres_sst
534         WRITE(numout,*) '      initial snow thickness in the north                    rn_hts_ini_n    = ', rn_hts_ini_n
535         WRITE(numout,*) '      initial snow thickness in the south                    rn_hts_ini_s    = ', rn_hts_ini_s 
536         WRITE(numout,*) '      initial ice thickness  in the north                    rn_hti_ini_n    = ', rn_hti_ini_n
537         WRITE(numout,*) '      initial ice thickness  in the south                    rn_hti_ini_s    = ', rn_hti_ini_s
538         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_n    = ', rn_ati_ini_n
539         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_s    = ', rn_ati_ini_s
540         WRITE(numout,*) '      initial  ice salinity  in the north                    rn_smi_ini_n    = ', rn_smi_ini_n
541         WRITE(numout,*) '      initial  ice salinity  in the south                    rn_smi_ini_s    = ', rn_smi_ini_s
542         WRITE(numout,*) '      initial  ice/snw temp  in the north                    rn_tmi_ini_n    = ', rn_tmi_ini_n
543         WRITE(numout,*) '      initial  ice/snw temp  in the south                    rn_tmi_ini_s    = ', rn_tmi_ini_s
544      ENDIF
545      !
546      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
547         !
548         ! set si structure
549         ALLOCATE( si(jpfldi), STAT=ierror )
550         IF( ierror > 0 ) THEN
551            CALL ctl_stop( 'Ice_ini in iceistate: unable to allocate si structure' )   ;   RETURN
552         ENDIF
553         !
554         DO ifpr = 1, jpfldi
555            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
556            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
557         END DO
558         !
559         ! fill si with slf_i and control print
560         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' )
561         !
562         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
563         !
564      ENDIF
565      !
566   END SUBROUTINE ice_istate_init
567
568#else
569   !!----------------------------------------------------------------------
570   !!   Default option :         Empty module         NO ESIM sea-ice model
571   !!----------------------------------------------------------------------
572#endif
573
574   !!======================================================================
575END MODULE iceistate
Note: See TracBrowser for help on using the repository browser.