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

Context Navigation

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

Last change on this file since 12725 was 12725, checked in by clem, 4 years ago
fix previous revision. Make sure a_ip and a_ip_eff are correct
Property svn:keywords set to `Id`
File size: 28.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
	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
	35	IMPLICIT NONE
	36	PRIVATE
	37
	38	PUBLIC ice_istate ! called by icestp.F90
	39	PUBLIC ice_istate_init ! called by icestp.F90
	40	!
[8813]	41	! !! namelist (namini)
[11536]	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
[12720]	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
[11536]	49	!
	50	! ! if ln_iceini_file = T
[12720]	51	INTEGER , PARAMETER :: jpfldi = 10 ! maximum number of files to read
[11536]	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)
[12720]	61	INTEGER , PARAMETER :: jp_hld = 10 ! index of pnd lid depth (m)
[11536]	62	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)
	63	!
[8586]	64	!!----------------------------------------------------------------------
[9598]	65	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	66	!! $Id$
[10413]	67	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[8586]	68	!!----------------------------------------------------------------------
	69	CONTAINS
	70
[11536]	71	SUBROUTINE ice_istate( kt )
[8586]	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
[11536]	90	!! where there is no ice
[8586]	91	!!--------------------------------------------------------------------
[11536]	92	INTEGER, INTENT(in) :: kt ! time step
	93	!!
[8813]	94	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[11536]	95	REAL(wp) :: ztmelts
[8586]	96	INTEGER , DIMENSION(4) :: itest
	97	REAL(wp), DIMENSION(jpi,jpj) :: z2d
	98	REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator
[11536]	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
[12720]	101	REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini, zhlid_ini !data from namelist or nc file
[11536]	102	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays
	103	!!
[12720]	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
[8586]	105	!--------------------------------------------------------------------
	106
	107	IF(lwp) WRITE(numout,*)
	108	IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
	109	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	110
[11536]	111	!---------------------------
	112	! 1) 1st init. of the fields
	113	!---------------------------
[8586]	114	!
[11536]	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
[8586]	120	DO jl = 1, jpl
[9416]	121	t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface
[10534]	122	cnd_ice(:,:,jl) = 0._wp ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
[8586]	123	END DO
[8813]	124	!
[11536]	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,:)
[8586]	138
[11536]	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
[12720]	158	v_il (:,:,:) = 0._wp
[11536]	159	a_ip_frac(:,:,:) = 0._wp
[12720]	160	a_ip_eff (:,:,:) = 0._wp
[11536]	161	h_ip (:,:,:) = 0._wp
[12720]	162	h_il (:,:,:) = 0._wp
[11536]	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	!------------------------------------------------------------------------
[8586]	171	IF( ln_iceini ) THEN
	172	! !---------------!
	173	IF( ln_iceini_file )THEN ! Read a file !
	174	! !---------------!
[11536]	175	WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp
	176	ELSEWHERE ; zswitch(:,:) = 0._wp
	177	END WHERE
[8586]	178	!
[11536]	179	CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
[8586]	180	!
[11536]	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 -- !
[12398]	187	! if fields do not exist then set them to the values present in the namelist (except for temperatures)
[11536]	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	!
[12398]	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)
[11536]	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)
[12398]	221	!
[12720]	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	!
[12398]	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)
[11536]	231	zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
[12720]	232	zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1)
[11536]	233	!
	234	! change the switch for the following
	235	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
	236	ELSEWHERE ; zswitch(:,:) = 0._wp
[8586]	237	END WHERE
[8813]	238	! !---------------!
[8586]	239	ELSE ! Read namelist !
	240	! !---------------!
[11536]	241	! no ice if (sst - Tfreez) >= thresold
[8813]	242	WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
	243	ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
[8586]	244	END WHERE
[8813]	245	!
[8586]	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(:,:)
[11536]	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(:,:)
[12720]	257	zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:)
[8586]	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(:,:)
[11536]	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(:,:)
[12720]	268	zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:)
[8586]	269	END WHERE
	270	!
	271	ENDIF
[11536]	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
[12720]	277	zhlid_ini(:,:) = 0._wp
[11536]	278	ENDIF
[8586]	279
[11536]	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
[8586]	293
[11536]	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 )
[12720]	304	CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti) , zhlid_ini )
[10413]	305
[11536]	306	! allocate temporary arrays
[12720]	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) )
[11536]	310
	311	! distribute 1-cat into jpl-cat: (jpijpj) -> (jpijpj,jpl)
[12720]	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 )
[8586]	318
[11536]	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 )
[12720]	333	CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d , h_il )
[8586]	334
[11536]	335	! deallocate temporary arrays
	336	DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
[12720]	337	& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d )
[11536]	338
	339	! calculate extensive and intensive variables
	340	CALL ice_var_salprof ! for sz_i
	341	DO jl = 1, jpl
[8586]	342	DO jj = 1, jpj
	343	DO ji = 1, jpi
[11536]	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)
[8586]	347	END DO
	348	END DO
	349	END DO
[8813]	350	!
[11536]	351	DO jl = 1, jpl
	352	DO jk = 1, nlay_s
[8586]	353	DO jj = 1, jpj
	354	DO ji = 1, jpi
[11536]	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 )
[8586]	358	END DO
	359	END DO
	360	END DO
	361	END DO
[8813]	362	!
[11536]	363	DO jl = 1, jpl
	364	DO jk = 1, nlay_i
[8586]	365	DO jj = 1, jpj
	366	DO ji = 1, jpi
[11536]	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 ) )
[8586]	373	END DO
	374	END DO
	375	END DO
	376	END DO
[8933]	377
[11536]	378	! Melt ponds
[12725]	379	WHERE( a_i > epsi10 ) ; a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
	380	ELSEWHERE ; a_ip_frac(:,:,:) = 0._wp
[11536]	381	END WHERE
[12720]	382	a_ip_eff(:,:,:) = a_ip_frac(:,:,:)
[11536]	383	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
[12720]	384	v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
[11536]	385
	386	! specific temperatures for coupled runs
	387	tn_ice(:,:,:) = t_su(:,:,:)
	388	t1_ice(:,:,:) = t_i (:,:,1,:)
[8813]	389	!
[8586]	390	ENDIF ! ln_iceini
[8813]	391	!
[11536]	392	at_i(:,:) = SUM( a_i, dim=3 )
[8586]	393	!
	394	!----------------------------------------------
[11536]	395	! 3) Snow-ice mass (case ice is fully embedded)
[8586]	396	!----------------------------------------------
[9935]	397	snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass
[8586]	398	snwice_mass_b(:,:) = snwice_mass(:,:)
	399	!
	400	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
[8813]	401	!
[8586]	402	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	403	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
[8813]	404	!
[8586]	405	IF( .NOT.ln_linssh ) THEN
[8813]	406	!
[8586]	407	WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
	408	ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE
[8813]	409	!
[8586]	410	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	411	e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
	412	e3t_b(:,:,jk) = e3t_n(:,:,jk)
	413	e3t_a(:,:,jk) = e3t_n(:,:,jk)
	414	END DO
[8813]	415	!
[8586]	416	! Reconstruction of all vertical scale factors at now and before time-steps
	417	! =========================================================================
	418	! Horizontal scale factor interpolations
	419	! --------------------------------------
	420	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
	421	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
	422	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	423	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
	424	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
	425	! Vertical scale factor interpolations
	426	! ------------------------------------
	427	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
	428	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	429	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
	430	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	431	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
	432	! t- and w- points depth
	433	! ----------------------
	434	!!gm not sure of that....
	435	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	436	gdepw_n(:,:,1) = 0.0_wp
	437	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
	438	DO jk = 2, jpk
	439	gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk )
	440	gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
	441	gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:)
	442	END DO
	443	ENDIF
	444	ENDIF
	445
	446	!------------------------------------
[11536]	447	! 4) store fields at before time-step
[8586]	448	!------------------------------------
	449	! it is only necessary for the 1st interpolation by Agrif
	450	a_i_b (:,:,:) = a_i (:,:,:)
	451	e_i_b (:,:,:,:) = e_i (:,:,:,:)
	452	v_i_b (:,:,:) = v_i (:,:,:)
	453	v_s_b (:,:,:) = v_s (:,:,:)
	454	e_s_b (:,:,:,:) = e_s (:,:,:,:)
	455	sv_i_b (:,:,:) = sv_i (:,:,:)
	456	u_ice_b(:,:) = u_ice(:,:)
	457	v_ice_b(:,:) = v_ice(:,:)
[10527]	458	! total concentration is needed for Lupkes parameterizations
	459	at_i_b (:,:) = at_i (:,:)
[8586]	460
[8885]	461	!!clem: output of initial state should be written here but it is impossible because
	462	!! the ocean and ice are in the same file
[10425]	463	!! CALL dia_wri_state( 'output.init' )
[8813]	464	!
[8586]	465	END SUBROUTINE ice_istate
	466
[8813]	467
[8586]	468	SUBROUTINE ice_istate_init
	469	!!-------------------------------------------------------------------
	470	!! * ROUTINE ice_istate_init *
	471	!!
[8813]	472	!! ** Purpose : Definition of initial state of the ice
[8586]	473	!!
[8813]	474	!! ** Method : Read the namini namelist and check the parameter
	475	!! values called at the first timestep (nit000)
[8586]	476	!!
[8813]	477	!! ** input : Namelist namini
[8586]	478	!!
	479	!!-----------------------------------------------------------------------------
[11229]	480	INTEGER :: ios ! Local integer output status for namelist read
[8813]	481	INTEGER :: ifpr, ierror
[8586]	482	!
	483	CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
[12720]	484	TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd, sn_hld
[8586]	485	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
	486	!
[11536]	487	NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
	488	& rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
	489	& rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
	490	& rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
[12720]	491	& rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, rn_hld_ini_n, rn_hld_ini_s, &
	492	& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, sn_hld, cn_dir
[8586]	493	!!-----------------------------------------------------------------------------
	494	!
	495	REWIND( numnam_ice_ref ) ! Namelist namini in reference namelist : Ice initial state
	496	READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
[11536]	497	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' )
[8586]	498	REWIND( numnam_ice_cfg ) ! Namelist namini in configuration namelist : Ice initial state
	499	READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
[11536]	500	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' )
[8586]	501	IF(lwm) WRITE ( numoni, namini )
[8813]	502	!
[8586]	503	slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
[11536]	504	slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi
	505	slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms
[12720]	506	slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd ; slf_i(jp_hld) = sn_hld
[8586]	507	!
	508	IF(lwp) THEN ! control print
	509	WRITE(numout,*)
	510	WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
	511	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	512	WRITE(numout,*) ' Namelist namini:'
[11536]	513	WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini
	514	WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file
	515	WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
	516	IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
	517	WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s
	518	WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s
	519	WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s
	520	WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s
	521	WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s
	522	WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s
	523	WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s
	524	WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s
	525	WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s
[12720]	526	WRITE(numout,*) ' initial pnd lid depth in the north-south rn_hld_ini = ', rn_hld_ini_n,rn_hld_ini_s
[11536]	527	ENDIF
[8586]	528	ENDIF
[8813]	529	!
[8586]	530	IF( ln_iceini_file ) THEN ! Ice initialization using input file
	531	!
	532	! set si structure
	533	ALLOCATE( si(jpfldi), STAT=ierror )
	534	IF( ierror > 0 ) THEN
[11536]	535	CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN
[8586]	536	ENDIF
[8813]	537	!
[8586]	538	DO ifpr = 1, jpfldi
	539	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
[11536]	540	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
[8586]	541	END DO
[8813]	542	!
[8586]	543	! fill si with slf_i and control print
[11536]	544	CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
[8813]	545	!
[8586]	546	ENDIF
[8813]	547	!
[11536]	548	IF( .NOT.ln_pnd ) THEN
	549	rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
	550	rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
[12720]	551	rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
	552	CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 & rn_hld_ini = 0 when no ponds' )
[11536]	553	ENDIF
	554	!
[8586]	555	END SUBROUTINE ice_istate_init
	556
	557	#else
	558	!!----------------------------------------------------------------------
[9570]	559	!! Default option : Empty module NO SI3 sea-ice model
[8586]	560	!!----------------------------------------------------------------------
	561	#endif
	562
	563	!!======================================================================
	564	END MODULE iceistate

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

Download in other formats: