New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

iceistate.F90 in NEMO/trunk/src/ICE – NEMO

Context Navigation

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

Last change on this file since 13601 was 13472, checked in by smasson, 4 years ago
trunk: commit changes from r4.0-HEAD from 13284 to 13449, see #2523
Property svn:keywords set to `Id`
File size: 30.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: