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 @ 12720

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

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