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/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE/iceistate.F90 @ 13463

Last change on this file since 13463 was 13463, checked in by andmirek, 4 years ago

Ticket #2195:update to trunk 13461

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