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/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE – NEMO

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90 @ 12726

Last change on this file since 12726 was 12726, checked in by clem, 4 years ago

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