source: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceistate.F90 @ 12636

Last change on this file since 12636 was 12636, checked in by dancopsey, 8 months ago

Add:

  • Melt pond lids
  • Melt pond maximum area and thickness
  • Melt pond vertical flushing
  • Area contributing to melt ponds depends on total ice fraction
File size: 27.1 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   !!            4.0  !  2018     (many people)       SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_istate       :  initialization of diagnostics ice variables
15   !!   ice_istate_init  :  initialization of ice state and namelist read
16   !!----------------------------------------------------------------------
17   USE phycst         ! physical constant
18   USE oce            ! dynamics and tracers variables
19   USE dom_oce        ! ocean domain
20   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd 
21   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
22   USE eosbn2         ! equation of state
23   USE domvvl         ! Variable volume
24   USE ice            ! sea-ice: variables
25   USE ice1D          ! sea-ice: thermodynamics variables
26   USE icetab         ! sea-ice: 1D <==> 2D transformation
27   USE icevar         ! sea-ice: operations
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   !                             !! ** namelist (namini) **
42   LOGICAL, PUBLIC  ::   ln_iceini        !: Ice initialization or not
43   LOGICAL, PUBLIC  ::   ln_iceini_file   !: Ice initialization from 2D netcdf file
44   REAL(wp) ::   rn_thres_sst
45   REAL(wp) ::   rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
46   REAL(wp) ::   rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
47   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n
48   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s
49   !
50   !                              ! if ln_iceini_file = T
51   INTEGER , PARAMETER ::   jpfldi = 9           ! maximum number of files to read
52   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m)
53   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m)
54   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction     (-)
55   INTEGER , PARAMETER ::   jp_smi = 4           ! index of ice salinity     (g/kg)
56   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temperature  (K)
57   INTEGER , PARAMETER ::   jp_tsu = 6           ! index of ice surface temp (K)
58   INTEGER , PARAMETER ::   jp_tms = 7           ! index of snw temperature  (K)
59   INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-)
60   INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
62   !   
63   !!----------------------------------------------------------------------
64   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
65   !! $Id$
66   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE ice_istate( kt )
71      !!-------------------------------------------------------------------
72      !!                    ***  ROUTINE ice_istate  ***
73      !!
74      !! ** Purpose :   defined the sea-ice initial state
75      !!
76      !! ** Method  :   This routine will put some ice where ocean
77      !!                is at the freezing point, then fill in ice
78      !!                state variables using prescribed initial
79      !!                values in the namelist           
80      !!
81      !! ** Steps   :   1) Set initial surface and basal temperatures
82      !!                2) Recompute or read sea ice state variables
83      !!                3) Fill in the ice thickness distribution using gaussian
84      !!                4) Fill in space-dependent arrays for state variables
85      !!                5) snow-ice mass computation
86      !!                6) store before fields
87      !!
88      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
89      !!              where there is no ice
90      !!--------------------------------------------------------------------
91      INTEGER, INTENT(in) ::   kt   ! time step
92      !!
93      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
94      REAL(wp) ::   ztmelts
95      INTEGER , DIMENSION(4)           ::   itest
96      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
97      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
98      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file
99      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
100      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file
101      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
102      !!
103      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d
104      !--------------------------------------------------------------------
105
106      IF(lwp) WRITE(numout,*)
107      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
108      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
109
110      !---------------------------
111      ! 1) 1st init. of the fields
112      !---------------------------
113      !
114      ! basal temperature (considered at freezing point)   [Kelvin]
115      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
116      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
117      !
118      ! surface temperature and conductivity
119      DO jl = 1, jpl
120         t_su   (:,:,jl) = rt0 * tmask(:,:,1)  ! temp at the surface
121         cnd_ice(:,:,jl) = 0._wp               ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
122      END DO
123      !
124      ! ice and snw temperatures
125      DO jl = 1, jpl
126         DO jk = 1, nlay_i
127            t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
128         END DO
129         DO jk = 1, nlay_s
130            t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
131         END DO
132      END DO
133      !
134      ! specific temperatures for coupled runs
135      tn_ice (:,:,:) = t_i (:,:,1,:)
136      t1_ice (:,:,:) = t_i (:,:,1,:)
137
138      ! heat contents
139      e_i (:,:,:,:) = 0._wp
140      e_s (:,:,:,:) = 0._wp
141     
142      ! general fields
143      a_i (:,:,:) = 0._wp
144      v_i (:,:,:) = 0._wp
145      v_s (:,:,:) = 0._wp
146      sv_i(:,:,:) = 0._wp
147      oa_i(:,:,:) = 0._wp
148      !
149      h_i (:,:,:) = 0._wp
150      h_s (:,:,:) = 0._wp
151      s_i (:,:,:) = 0._wp
152      o_i (:,:,:) = 0._wp
153      !
154      ! melt ponds
155      a_ip     (:,:,:) = 0._wp
156      v_ip     (:,:,:) = 0._wp
157      lh_ip    (:,:,:) = 0._wp
158      a_ip_frac(:,:,:) = 0._wp
159      a_ip_eff (:,:,:) = 0._wp
160      h_ip     (:,:,:) = 0._wp
161      !
162      ! ice velocities
163      u_ice (:,:) = 0._wp
164      v_ice (:,:) = 0._wp
165      !
166      !------------------------------------------------------------------------
167      ! 2) overwrite some of the fields with namelist parameters or netcdf file
168      !------------------------------------------------------------------------
169      IF( ln_iceini ) THEN
170         !                             !---------------!
171         IF( ln_iceini_file )THEN      ! Read a file   !
172            !                          !---------------!
173            WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
174            ELSEWHERE                     ;   zswitch(:,:) = 0._wp
175            END WHERE
176            !
177            CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
178            !
179            ! -- mandatory fields -- !
180            zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
181            zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
182            zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
183
184            ! -- optional fields -- !
185            !    if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
186            !
187            ! ice salinity
188            IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
189               &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
190            zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
191            !
192            ! ice temperature
193            IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) &
194               &     si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
195            ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
196            !
197            ! surface temperature => set to ice temperature if it exists
198            IF    ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
199                     si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
200            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
201                     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
202            ENDIF
203            zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
204            !
205            ! snow temperature => set to ice temperature if it exists
206            IF    ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
207                     si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
208            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
209                     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
210            ENDIF
211            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1)
212            !
213            ! pond concentration
214            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
215               &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
216               &                              * si(jp_ati)%fnow(:,:,1) 
217            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
218            !
219            ! pond depth
220            IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
221               &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
222            zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
223            !
224            ! change the switch for the following
225            WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
226            ELSEWHERE                         ;   zswitch(:,:) = 0._wp
227            END WHERE
228            !                          !---------------!
229         ELSE                          ! Read namelist !
230            !                          !---------------!
231            ! no ice if (sst - Tfreez) >= thresold
232            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
233            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
234            END WHERE
235            !
236            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
237            WHERE( ff_t(:,:) >= 0._wp )
238               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
239               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
240               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
241               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
242               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
243               zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
244               ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
245               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
246               zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
247            ELSEWHERE
248               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
249               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
250               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
251               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
252               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
253               zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
254               ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
255               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
256               zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
257            END WHERE
258            !
259         ENDIF
260
261         ! make sure ponds = 0 if no ponds scheme
262         IF ( .NOT.ln_pnd ) THEN
263            zapnd_ini(:,:) = 0._wp
264            zhpnd_ini(:,:) = 0._wp
265         ENDIF
266         
267         !-------------!
268         ! fill fields !
269         !-------------!
270         ! select ice covered grid points
271         npti = 0 ; nptidx(:) = 0
272         DO jj = 1, jpj
273            DO ji = 1, jpi
274               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
275                  npti         = npti  + 1
276                  nptidx(npti) = (jj - 1) * jpi + ji
277               ENDIF
278            END DO
279         END DO
280
281         ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
282         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini )
283         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini )
284         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini )
285         CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
286         CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
287         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
288         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
289         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
290         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
291
292         ! allocate temporary arrays
293         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
294            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
295         
296         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
297         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   &
298            &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   &
299            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), &
300            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d )
301
302         ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
303         DO jl = 1, jpl
304            zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
305            zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
306         END DO
307         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    )
308         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    )
309         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    )
310         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d )
311         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d )
312         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   )
313         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    )
314         CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
315         CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
316
317         ! deallocate temporary arrays
318         DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
319            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
320
321         ! calculate extensive and intensive variables
322         CALL ice_var_salprof ! for sz_i
323         DO jl = 1, jpl
324            DO jj = 1, jpj
325               DO ji = 1, jpi
326                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
327                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
328                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
329               END DO
330            END DO
331         END DO
332         !
333         DO jl = 1, jpl
334            DO jk = 1, nlay_s
335               DO jj = 1, jpj
336                  DO ji = 1, jpi
337                     t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
338                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
339                        &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
340                  END DO
341               END DO
342            END DO
343         END DO
344         !
345         DO jl = 1, jpl
346            DO jk = 1, nlay_i
347               DO jj = 1, jpj
348                  DO ji = 1, jpi
349                     t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
350                     ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
351                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
352                        &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
353                        &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
354                        &                       - rcp   * ( ztmelts - rt0 ) )
355                  END DO
356               END DO
357            END DO
358         END DO
359
360         ! Melt ponds
361         WHERE( a_i > epsi10 )
362            a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
363         ELSEWHERE
364            a_ip_frac(:,:,:) = 0._wp
365         END WHERE
366         a_ip_eff(:,:,:) = a_ip_frac(:,:,:)
367         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
368         
369         ! specific temperatures for coupled runs
370         tn_ice(:,:,:) = t_su(:,:,:)
371         t1_ice(:,:,:) = t_i (:,:,1,:)
372         !
373      ENDIF ! ln_iceini
374      !
375      at_i(:,:) = SUM( a_i, dim=3 )
376      !
377      !----------------------------------------------
378      ! 3) Snow-ice mass (case ice is fully embedded)
379      !----------------------------------------------
380      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass
381      snwice_mass_b(:,:) = snwice_mass(:,:)
382      !
383      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
384         !
385         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
386         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
387         !
388         IF( .NOT.ln_linssh ) THEN
389            !
390            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
391            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
392            !
393            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
394               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
395               e3t_b(:,:,jk) = e3t_n(:,:,jk)
396               e3t_a(:,:,jk) = e3t_n(:,:,jk)
397            END DO
398            !
399            ! Reconstruction of all vertical scale factors at now and before time-steps
400            ! =========================================================================
401            ! Horizontal scale factor interpolations
402            ! --------------------------------------
403            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
404            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
405            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
406            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
407            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
408            ! Vertical scale factor interpolations
409            ! ------------------------------------
410            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
411            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
412            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
413            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
414            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
415            ! t- and w- points depth
416            ! ----------------------
417            !!gm not sure of that....
418            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
419            gdepw_n(:,:,1) = 0.0_wp
420            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
421            DO jk = 2, jpk
422               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
423               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
424               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
425            END DO
426         ENDIF
427      ENDIF
428     
429      !------------------------------------
430      ! 4) store fields at before time-step
431      !------------------------------------
432      ! it is only necessary for the 1st interpolation by Agrif
433      a_i_b  (:,:,:)   = a_i  (:,:,:)
434      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
435      v_i_b  (:,:,:)   = v_i  (:,:,:)
436      v_s_b  (:,:,:)   = v_s  (:,:,:)
437      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
438      sv_i_b (:,:,:)   = sv_i (:,:,:)
439      oa_i_b (:,:,:)   = oa_i (:,:,:)
440      u_ice_b(:,:)     = u_ice(:,:)
441      v_ice_b(:,:)     = v_ice(:,:)
442      ! total concentration is needed for Lupkes parameterizations
443      at_i_b (:,:)     = at_i (:,:) 
444
445!!clem: output of initial state should be written here but it is impossible because
446!!      the ocean and ice are in the same file
447!!      CALL dia_wri_state( 'output.init' )
448      !
449   END SUBROUTINE ice_istate
450
451
452   SUBROUTINE ice_istate_init
453      !!-------------------------------------------------------------------
454      !!                   ***  ROUTINE ice_istate_init  ***
455      !!       
456      !! ** Purpose :   Definition of initial state of the ice
457      !!
458      !! ** Method  :   Read the namini namelist and check the parameter
459      !!              values called at the first timestep (nit000)
460      !!
461      !! ** input   :  Namelist namini
462      !!
463      !!-----------------------------------------------------------------------------
464      INTEGER ::   ios   ! Local integer output status for namelist read
465      INTEGER ::   ifpr, ierror
466      !
467      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
468      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
469      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
470      !
471      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
472         &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
473         &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
474         &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
475         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
476         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
477      !!-----------------------------------------------------------------------------
478      !
479      REWIND( numnam_ice_ref )              ! Namelist namini in reference namelist : Ice initial state
480      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
481901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namini in reference namelist' )
482      REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state
483      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
484902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namini in configuration namelist' )
485      IF(lwm) WRITE ( numoni, namini )
486      !
487      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
488      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi
489      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms
490      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd
491      !
492      IF(lwp) THEN                          ! control print
493         WRITE(numout,*)
494         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
495         WRITE(numout,*) '~~~~~~~~~~~~~~~'
496         WRITE(numout,*) '   Namelist namini:'
497         WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini
498         WRITE(numout,*) '      ice initialization from a netcdf file            ln_iceini_file = ', ln_iceini_file
499         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst
500         IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
501            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
502            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s
503            WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s
504            WRITE(numout,*) '      initial ice salinity  in the north-south         rn_smi_ini     = ', rn_smi_ini_n,rn_smi_ini_s
505            WRITE(numout,*) '      initial surf temperat in the north-south         rn_tsu_ini     = ', rn_tsu_ini_n,rn_tsu_ini_s
506            WRITE(numout,*) '      initial ice temperat  in the north-south         rn_tmi_ini     = ', rn_tmi_ini_n,rn_tmi_ini_s
507            WRITE(numout,*) '      initial snw temperat  in the north-south         rn_tms_ini     = ', rn_tms_ini_n,rn_tms_ini_s
508            WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s
509            WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s
510         ENDIF
511      ENDIF
512      !
513      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
514         !
515         ! set si structure
516         ALLOCATE( si(jpfldi), STAT=ierror )
517         IF( ierror > 0 ) THEN
518            CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' )   ;   RETURN
519         ENDIF
520         !
521         DO ifpr = 1, jpfldi
522            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
523            IF( slf_i(ifpr)%ln_tint )  ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
524         END DO
525         !
526         ! fill si with slf_i and control print
527         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
528         !
529      ENDIF
530      !
531      IF( .NOT.ln_pnd ) THEN
532         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
533         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
534         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
535      ENDIF
536      !
537   END SUBROUTINE ice_istate_init
538
539#else
540   !!----------------------------------------------------------------------
541   !!   Default option :         Empty module         NO SI3 sea-ice model
542   !!----------------------------------------------------------------------
543#endif
544
545   !!======================================================================
546END MODULE iceistate
Note: See TracBrowser for help on using the repository browser.