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 NEMO/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/iceistate.F90 @ 13418

Last change on this file since 13418 was 13418, checked in by techene, 4 years ago

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