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

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

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

Download in other formats: