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.

sbcice_cice.F90 in branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9695

Last change on this file since 9695 was 9695, checked in by dancopsey, 6 years ago
Apply -Lsub to all cases of zevap_ice. Multiply v_ip by a_i.
File size: 46.3 KB

Rev	Line
[2874]	1	MODULE sbcice_cice
	2	!!======================================================================
	3	!! * MODULE sbcice_cice *
	4	!! To couple with sea ice model CICE (LANL)
	5	!!=====================================================================
	6	#if defined key_cice
	7	!!----------------------------------------------------------------------
	8	!! 'key_cice' : CICE sea-ice model
	9	!!----------------------------------------------------------------------
	10	!! sbc_ice_cice : sea-ice model time-stepping and update ocean sbc over ice-covered area
	11	!!----------------------------------------------------------------------
	12	USE oce ! ocean dynamics and tracers
	13	USE dom_oce ! ocean space and time domain
[3275]	14	USE domvvl
[9500]	15	USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater
	16	USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0
[2874]	17	USE in_out_manager ! I/O manager
[4990]	18	USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit
[2874]	19	USE lib_mpp ! distributed memory computing library
	20	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	21	USE daymod ! calendar
	22	USE fldread ! read input fields
	23	USE sbc_oce ! Surface boundary condition: ocean fields
	24	USE sbc_ice ! Surface boundary condition: ice fields
[7646]	25	USE sbcblk ! Surface boundary condition: bulk
[2874]	26	USE sbccpl
	27
	28	USE ice_kinds_mod
	29	USE ice_blocks
	30	USE ice_domain
	31	USE ice_domain_size
	32	USE ice_boundary
	33	USE ice_constants
	34	USE ice_gather_scatter
	35	USE ice_calendar, only: dt
[9500]	36	# if defined key_cice4
[9498]	37	USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
[2874]	38	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[5133]	39	strocnxT,strocnyT, &
[3189]	40	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, &
	41	fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, &
[2874]	42	flatn_f,fsurfn_f,fcondtopn_f, &
	43	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
[9500]	44	swvdr,swvdf,swidr,swidf,Tf
[4990]	45	USE ice_therm_vertical, only: calc_Tsfc
	46	#else
[9500]	47	USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,&
	48	vsnon,vice,vicen,nt_Tsfc
[4990]	49	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[5133]	50	strocnxT,strocnyT, &
[9500]	51	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai, &
	52	fresh_ai,fhocn_ai,fswthru_ai,frzmlt, &
[4990]	53	flatn_f,fsurfn_f,fcondtopn_f, &
[9500]	54	#ifdef key_asminc
	55	daice_da,fresh_da,fsalt_da, &
	56	#endif
[4990]	57	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
[9500]	58	swvdr,swvdf,swidr,swidf,Tf, &
	59	!! When using NEMO with CICE, this change requires use of
	60	!! one of the following two CICE branches:
	61	!! - at CICE5.0, hadax/r1015_GSI8_with_GSI7
	62	!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
	63	keffn_top,Tn_top
	64
	65	USE ice_therm_shared, only: calc_Tsfc, heat_capacity
	66	USE ice_shortwave, only: apeffn
[4990]	67	#endif
[2874]	68	USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
[3176]	69	USE ice_atmo, only: calc_strair
[2874]	70
	71	USE CICE_InitMod
	72	USE CICE_RunMod
	73	USE CICE_FinalMod
	74
	75	IMPLICIT NONE
	76	PRIVATE
	77
	78	PUBLIC cice_sbc_init ! routine called by sbc_init
	79	PUBLIC cice_sbc_final ! routine called by sbc_final
	80	PUBLIC sbc_ice_cice ! routine called by sbc
	81
[4627]	82	INTEGER :: ji_off
	83	INTEGER :: jj_off
[3625]	84
[2874]	85	INTEGER , PARAMETER :: jpfld = 13 ! maximum number of files to read
	86	INTEGER , PARAMETER :: jp_snow = 1 ! index of snow file
	87	INTEGER , PARAMETER :: jp_rain = 2 ! index of rain file
	88	INTEGER , PARAMETER :: jp_sblm = 3 ! index of sublimation file
	89	INTEGER , PARAMETER :: jp_top1 = 4 ! index of category 1 topmelt file
	90	INTEGER , PARAMETER :: jp_top2 = 5 ! index of category 2 topmelt file
	91	INTEGER , PARAMETER :: jp_top3 = 6 ! index of category 3 topmelt file
	92	INTEGER , PARAMETER :: jp_top4 = 7 ! index of category 4 topmelt file
	93	INTEGER , PARAMETER :: jp_top5 = 8 ! index of category 5 topmelt file
	94	INTEGER , PARAMETER :: jp_bot1 = 9 ! index of category 1 botmelt file
	95	INTEGER , PARAMETER :: jp_bot2 = 10 ! index of category 2 botmelt file
	96	INTEGER , PARAMETER :: jp_bot3 = 11 ! index of category 3 botmelt file
	97	INTEGER , PARAMETER :: jp_bot4 = 12 ! index of category 4 botmelt file
	98	INTEGER , PARAMETER :: jp_bot5 = 13 ! index of category 5 botmelt file
	99	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
	100
	101	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE :: png ! local array used in sbc_cice_ice
	102
[5836]	103	!!----------------------------------------------------------------------
	104	!! NEMO/OPA 3.7 , NEMO-consortium (2015)
[5215]	105	!! $Id$
[5836]	106	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	107	!!----------------------------------------------------------------------
[2874]	108	CONTAINS
	109
	110	INTEGER FUNCTION sbc_ice_cice_alloc()
	111	!!----------------------------------------------------------------------
	112	!! * FUNCTION sbc_ice_cice_alloc *
	113	!!----------------------------------------------------------------------
	114	ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
	115	IF( lk_mpp ) CALL mpp_sum ( sbc_ice_cice_alloc )
	116	IF( sbc_ice_cice_alloc > 0 ) CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
	117	END FUNCTION sbc_ice_cice_alloc
	118
[4990]	119	SUBROUTINE sbc_ice_cice( kt, ksbc )
[2874]	120	!!---------------------------------------------------------------------
	121	!! * ROUTINE sbc_ice_cice *
	122	!!
	123	!! ** Purpose : update the ocean surface boundary condition via the
	124	!! CICE Sea Ice Model time stepping
	125	!!
[3040]	126	!! ** Method : - Get any extra forcing fields for CICE
	127	!! - Prepare forcing fields
[2874]	128	!! - CICE model time stepping
	129	!! - call the routine that computes mass and
	130	!! heat fluxes at the ice/ocean interface
	131	!!
	132	!! ** Action : - time evolution of the CICE sea-ice model
	133	!! - update all sbc variables below sea-ice:
[3625]	134	!! utau, vtau, qns , qsr, emp , sfx
[2874]	135	!!---------------------------------------------------------------------
	136	INTEGER, INTENT(in) :: kt ! ocean time step
[4990]	137	INTEGER, INTENT(in) :: ksbc ! surface forcing type
[2874]	138	!!----------------------------------------------------------------------
[3193]	139	!
[2874]	140	! !----------------------!
	141	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Ice time-step only !
	142	! !----------------------!
	143
	144	! Make sure any fluxes required for CICE are set
[4990]	145	IF ( ksbc == jp_flx ) THEN
[2874]	146	CALL cice_sbc_force(kt)
[5407]	147	ELSE IF ( ksbc == jp_purecpl ) THEN
[9019]	148	CALL sbc_cpl_ice_flx( fr_i )
[2874]	149	ENDIF
	150
[4990]	151	CALL cice_sbc_in ( kt, ksbc )
[2874]	152	CALL CICE_Run
[4990]	153	CALL cice_sbc_out ( kt, ksbc )
[2874]	154
[5407]	155	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(kt+1)
[2874]	156
	157	ENDIF ! End sea-ice time step only
[3193]	158	!
[2874]	159	END SUBROUTINE sbc_ice_cice
	160
[5836]	161
	162	SUBROUTINE cice_sbc_init( ksbc )
[2874]	163	!!---------------------------------------------------------------------
	164	!! * ROUTINE cice_sbc_init *
[3040]	165	!! ** Purpose: Initialise ice related fields for NEMO and coupling
[2874]	166	!!
[5836]	167	!!---------------------------------------------------------------------
[4990]	168	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
[9125]	169	REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
[9500]	170	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d
[4990]	171	INTEGER :: ji, jj, jl, jk ! dummy loop indices
[2874]	172	!!---------------------------------------------------------------------
[3193]	173	!
[2874]	174	IF(lwp) WRITE(numout,*)'cice_sbc_init'
	175
[4627]	176	ji_off = INT ( (jpiglo - nx_global) / 2 )
	177	jj_off = INT ( (jpjglo - ny_global) / 2 )
	178
[9500]	179	! Initialize CICE
[3176]	180	CALL CICE_Initialize
[2874]	181
[9500]	182	! Do some CICE consistency checks
[5407]	183	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
[3193]	184	IF ( calc_strair .OR. calc_Tsfc ) THEN
	185	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
	186	ENDIF
[7646]	187	ELSEIF (ksbc == jp_blk) THEN
[3193]	188	IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
	189	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
	190	ENDIF
	191	ENDIF
[3176]	192
	193
[9500]	194	! allocate sbc_ice and sbc_cice arrays
	195	IF( sbc_ice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_alloc : unable to allocate arrays' )
[2874]	196	IF( sbc_ice_cice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
	197
[9500]	198	! Ensure that no temperature points are below freezing if not a NEMO restart
[2874]	199	IF( .NOT. ln_rstart ) THEN
[9500]	200
	201	DO jk=1,jpk
	202	CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk) )
	203	ENDDO
	204	tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d )
[2874]	205	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
	206
[9500]	207	#if defined key_nemocice_decomp
	208	! Pass initial SST from NEMO to CICE so ice is initialised correctly if
	209	! there is no restart file.
	210	! Values from a CICE restart file would overwrite this
	211	CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)
	212	#endif
[2874]	213
[9500]	214	ENDIF
	215
	216	! calculate surface freezing temperature and send to CICE
	217	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1))
	218	CALL nemo2cice(sstfrz,Tf, 'T', 1. )
	219
[3193]	220	CALL cice2nemo(aice,fr_i, 'T', 1. )
[5407]	221	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
[3625]	222	DO jl=1,ncat
	223	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	224	ENDDO
	225	ENDIF
[2874]	226
	227	! T point to U point
	228	! T point to V point
[9500]	229	fr_iu(:,:)=0.0
	230	fr_iv(:,:)=0.0
[3193]	231	DO jj=1,jpjm1
	232	DO ji=1,jpim1
	233	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	234	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	235	ENDDO
	236	ENDDO
[2874]	237
[9098]	238	CALL lbc_lnk_multi( fr_iu , 'U', 1., fr_iv , 'V', 1. )
[3625]	239
[9019]	240	! set the snow+ice mass
	241	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	242	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	243	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	244	snwice_mass_b(:,:) = snwice_mass(:,:)
	245
[6140]	246	IF( .NOT.ln_rstart ) THEN
[9019]	247	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
[4990]	248	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	249	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
[6140]	250
	251	!!gm This should be put elsewhere.... (same remark for limsbc)
	252	!!gm especially here it is assumed zstar coordinate, but it can be ztilde....
	253	IF( .NOT.ln_linssh ) THEN
	254	!
	255	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	256	e3t_n(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshn(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
	257	e3t_b(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshb(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
	258	ENDDO
	259	e3t_a(:,:,:) = e3t_b(:,:,:)
	260	! Reconstruction of all vertical scale factors at now and before time-steps
	261	! =============================================================================
	262	! Horizontal scale factor interpolations
	263	! --------------------------------------
	264	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
	265	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
	266	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	267	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
	268	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
	269	! Vertical scale factor interpolations
	270	! ------------------------------------
	271	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
	272	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	273	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
	274	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	275	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
	276	! t- and w- points depth
	277	! ----------------------
	278	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	279	gdepw_n(:,:,1) = 0.0_wp
	280	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
	281	DO jk = 2, jpk
	282	gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk)
	283	gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
	284	gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:)
	285	END DO
	286	ENDIF
[4990]	287	ENDIF
[3625]	288	ENDIF
[9500]	289
	290	#if defined key_asminc
	291	! Initialize fresh water and salt fluxes from data assim
	292	! and data assimilation index to cice
	293	nfresh_da(:,:) = 0.0
	294	nfsalt_da(:,:) = 0.0
	295	ndaice_da(:,:) = 0.0
	296	#endif
[6140]	297	!
[9500]	298	! In coupled mode get extra fields from CICE for passing back to atmosphere
	299
	300	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(nit000)
	301	!
[2874]	302	END SUBROUTINE cice_sbc_init
	303
[3152]	304
[5836]	305	SUBROUTINE cice_sbc_in( kt, ksbc )
[2874]	306	!!---------------------------------------------------------------------
	307	!! * ROUTINE cice_sbc_in *
[3040]	308	!! ** Purpose: Set coupling fields and pass to CICE
[2874]	309	!!---------------------------------------------------------------------
[3152]	310	INTEGER, INTENT(in ) :: kt ! ocean time step
[4990]	311	INTEGER, INTENT(in ) :: ksbc ! surface forcing type
[5836]	312	!
[3625]	313	INTEGER :: ji, jj, jl ! dummy loop indices
[9125]	314	REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zpice
	315	REAL(wp), DIMENSION(jpi,jpj,ncat) :: ztmpn
[3625]	316	REAL(wp) :: zintb, zintn ! dummy argument
[3152]	317	!!---------------------------------------------------------------------
[3193]	318	!
	319	IF( kt == nit000 ) THEN
[2874]	320	IF(lwp) WRITE(numout,*)'cice_sbc_in'
[3193]	321	ENDIF
[2874]	322
[3193]	323	ztmp(:,:)=0.0
[2874]	324
	325	! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
	326	! the first time-step)
	327
	328	! forced and coupled case
	329
[5407]	330	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
[2874]	331
[3193]	332	ztmpn(:,:,:)=0.0
[2874]	333
	334	! x comp of wind stress (CI_1)
	335	! U point to F point
[3193]	336	DO jj=1,jpjm1
	337	DO ji=1,jpi
	338	ztmp(ji,jj) = 0.5 * ( fr_iu(ji,jj) * utau(ji,jj) &
	339	+ fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
	340	ENDDO
	341	ENDDO
	342	CALL nemo2cice(ztmp,strax,'F', -1. )
[2874]	343
	344	! y comp of wind stress (CI_2)
	345	! V point to F point
[3193]	346	DO jj=1,jpj
	347	DO ji=1,jpim1
	348	ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
	349	+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
	350	ENDDO
	351	ENDDO
	352	CALL nemo2cice(ztmp,stray,'F', -1. )
[2874]	353
[9500]	354
	355	! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in
	356	! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby
	357	! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing
	358	! gridbox mean fluxes in the UM by future ice concentration obtained through
	359	! OASIS. This allows for a much more realistic apportionment of energy through
	360	! the ice - and conserves energy.
	361	! Therefore the fluxes are now divided by ice concentration in the coupled
	362	! formulation (jp_purecpl) as well as for jp_flx. This NEMO branch should only
	363	! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at
	364	! which point the GSI8 UM changes were committed.
	365
[2874]	366	! Surface downward latent heat flux (CI_5)
[4990]	367	IF (ksbc == jp_flx) THEN
[3625]	368	DO jl=1,ncat
[9695]	369	ztmpn(:,:,jl)=zevap_ice(:,:,1)a_i(:,:,jl) (-1.0*Lsub)
[3193]	370	ENDDO
[9500]	371	ELSE IF (ksbc == jp_purecpl) THEN
	372	DO jl=1,ncat
[9675]	373	ztmpn(:,:,jl)=zevap_ice(:,:,jl)a_i(:,:,jl) (-1.0*Lsub)
[3193]	374	ENDDO
[9500]	375	ELSE
[9675]	376	!In coupled mode - zevap_ice calculated in sbc_cpl for each category
[9695]	377	ztmpn(:,:,1:ncat)=zevap_ice(:,:,1:ncat) * (-1.0*Lsub)
[3193]	378	ENDIF
[9500]	379
[3625]	380	DO jl=1,ncat
	381	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
[2874]	382
	383	! GBM conductive flux through ice (CI_6)
	384	! Convert to GBM
[9500]	385	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
[9675]	386	ztmp(:,:) = qcn_ice(:,:,jl)*a_i(:,:,jl)
[3193]	387	ELSE
[9675]	388	ztmp(:,:) = qcn_ice(:,:,jl)
[3193]	389	ENDIF
[3625]	390	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
[2874]	391
	392	! GBM surface heat flux (CI_7)
	393	! Convert to GBM
[9500]	394	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
[9675]	395	ztmp(:,:) = (qml_ice(:,:,jl)+qcn_ice(:,:,jl))*a_i(:,:,jl)
[3193]	396	ELSE
[9675]	397	ztmp(:,:) = (qml_ice(:,:,jl)+qcn_ice(:,:,jl))
[3193]	398	ENDIF
[3625]	399	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
[3193]	400	ENDDO
[2874]	401
[7646]	402	ELSE IF (ksbc == jp_blk) THEN
[2874]	403
[7646]	404	! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
[2874]	405	! x comp and y comp of atmosphere surface wind (CICE expects on T points)
[3193]	406	ztmp(:,:) = wndi_ice(:,:)
	407	CALL nemo2cice(ztmp,uatm,'T', -1. )
	408	ztmp(:,:) = wndj_ice(:,:)
	409	CALL nemo2cice(ztmp,vatm,'T', -1. )
	410	ztmp(:,:) = SQRT ( wndi_ice(:,:)2 + wndj_ice(:,:)2 )
	411	CALL nemo2cice(ztmp,wind,'T', 1. ) ! Wind speed (m/s)
	412	ztmp(:,:) = qsr_ice(:,:,1)
	413	CALL nemo2cice(ztmp,fsw,'T', 1. ) ! Incoming short-wave (W/m^2)
	414	ztmp(:,:) = qlw_ice(:,:,1)
	415	CALL nemo2cice(ztmp,flw,'T', 1. ) ! Incoming long-wave (W/m^2)
	416	ztmp(:,:) = tatm_ice(:,:)
	417	CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K)
	418	CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K)
[2874]	419	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
[3193]	420	ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
	421	! Constant (101000.) atm pressure assumed
	422	CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3)
	423	ztmp(:,:) = qatm_ice(:,:)
	424	CALL nemo2cice(ztmp,Qa,'T', 1. ) ! Specific humidity (kg/kg)
	425	ztmp(:,:)=10.0
	426	CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m)
[2874]	427
	428	! May want to check all values are physically realistic (as in CICE routine
	429	! prepare_forcing)?
	430
	431	! Divide shortwave into spectral bands (as in prepare_forcing)
[3193]	432	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
[2874]	433	CALL nemo2cice(ztmp,swvdr,'T', 1. )
[3193]	434	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
[2874]	435	CALL nemo2cice(ztmp,swvdf,'T', 1. )
[3193]	436	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
[2874]	437	CALL nemo2cice(ztmp,swidr,'T', 1. )
[3193]	438	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
[2874]	439	CALL nemo2cice(ztmp,swidf,'T', 1. )
	440
	441	ENDIF
	442
[9500]	443	#if defined key_asminc
	444	!Ice concentration change (from assimilation)
	445	ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
	446	Call nemo2cice(ztmp,daice_da,'T', 1. )
	447	#endif
	448
[2874]	449	! Snowfall
[4990]	450	! Ensure fsnow is positive (as in CICE routine prepare_forcing)
	451	IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
[3193]	452	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
	453	CALL nemo2cice(ztmp,fsnow,'T', 1. )
[2874]	454
	455	! Rainfall
[4990]	456	IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
[3193]	457	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
	458	CALL nemo2cice(ztmp,frain,'T', 1. )
[2874]	459
[9500]	460	! Recalculate freezing temperature and send to CICE
	461	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1))
	462	CALL nemo2cice(sstfrz,Tf,'T', 1. )
	463
[2874]	464	! Freezing/melting potential
[3275]	465	! Calculated over NEMO leapfrog timestep (hence 2*dt)
[9500]	466	nfrzmlt(:,:)=rau0rcpe3t_m(:,:)(sstfrz(:,:)-sst_m(:,:))/(2.0dt)
	467	CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. )
[2874]	468
	469	! SST and SSS
	470
[3193]	471	CALL nemo2cice(sst_m,sst,'T', 1. )
	472	CALL nemo2cice(sss_m,sss,'T', 1. )
[2874]	473
[9500]	474	IF( ksbc == jp_purecpl ) THEN
	475	! Sea ice surface skin temperature
	476	DO jl=1,ncat
	477	CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.)
	478	ENDDO
	479	ENDIF
	480
[2874]	481	! x comp and y comp of surface ocean current
	482	! U point to F point
[3193]	483	DO jj=1,jpjm1
	484	DO ji=1,jpi
	485	ztmp(ji,jj)=0.5(ssu_m(ji,jj)+ssu_m(ji,jj+1))fmask(ji,jj,1)
	486	ENDDO
	487	ENDDO
	488	CALL nemo2cice(ztmp,uocn,'F', -1. )
[2874]	489
	490	! V point to F point
[3193]	491	DO jj=1,jpj
	492	DO ji=1,jpim1
	493	ztmp(ji,jj)=0.5(ssv_m(ji,jj)+ssv_m(ji+1,jj))fmask(ji,jj,1)
	494	ENDDO
	495	ENDDO
	496	CALL nemo2cice(ztmp,vocn,'F', -1. )
[2874]	497
[9019]	498	IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==!
[3625]	499	!
	500	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	501	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
	502	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
	503	!
	504	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	505	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
	506	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
	507	!
	508	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
	509	!
	510	!
	511	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	512	zpice(:,:) = ssh_m(:,:)
	513	ENDIF
	514
[3189]	515	! x comp and y comp of sea surface slope (on F points)
	516	! T point to F point
[5836]	517	DO jj = 1, jpjm1
	518	DO ji = 1, jpim1
	519	ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj )) * r1_e1u(ji,jj ) &
	520	& + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1) ) * fmask(ji,jj,1)
	521	END DO
	522	END DO
	523	CALL nemo2cice( ztmp,ss_tltx,'F', -1. )
[3189]	524
	525	! T point to F point
[5836]	526	DO jj = 1, jpjm1
	527	DO ji = 1, jpim1
	528	ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj)) * r1_e2v(ji ,jj) &
	529	& + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj) ) * fmask(ji,jj,1)
	530	END DO
	531	END DO
[3193]	532	CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
	533	!
[2874]	534	END SUBROUTINE cice_sbc_in
	535
[3152]	536
[5836]	537	SUBROUTINE cice_sbc_out( kt, ksbc )
[2874]	538	!!---------------------------------------------------------------------
	539	!! * ROUTINE cice_sbc_out *
[3040]	540	!! ** Purpose: Get fields from CICE and set surface fields for NEMO
[3152]	541	!!---------------------------------------------------------------------
[2874]	542	INTEGER, INTENT( in ) :: kt ! ocean time step
[4990]	543	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
[3152]	544
[3625]	545	INTEGER :: ji, jj, jl ! dummy loop indices
[9125]	546	REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
[2874]	547	!!---------------------------------------------------------------------
[3193]	548	!
[3152]	549	IF( kt == nit000 ) THEN
[2874]	550	IF(lwp) WRITE(numout,*)'cice_sbc_out'
[3152]	551	ENDIF
	552
[2874]	553	! x comp of ocean-ice stress
[3625]	554	CALL cice2nemo(strocnx,ztmp1,'F', -1. )
[3193]	555	ss_iou(:,:)=0.0
[2874]	556	! F point to U point
[3193]	557	DO jj=2,jpjm1
	558	DO ji=2,jpim1
[3625]	559	ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
[3193]	560	ENDDO
	561	ENDDO
	562	CALL lbc_lnk( ss_iou , 'U', -1. )
[2874]	563
	564	! y comp of ocean-ice stress
[3625]	565	CALL cice2nemo(strocny,ztmp1,'F', -1. )
[3193]	566	ss_iov(:,:)=0.0
[2874]	567	! F point to V point
	568
[3193]	569	DO jj=1,jpjm1
	570	DO ji=2,jpim1
[3625]	571	ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
[3193]	572	ENDDO
	573	ENDDO
	574	CALL lbc_lnk( ss_iov , 'V', -1. )
[2874]	575
	576	! x and y comps of surface stress
	577	! Combine wind stress and ocean-ice stress
	578	! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
[5133]	579	! strocnx and strocny already weighted by ice fraction in CICE so not done here
[2874]	580
[3193]	581	utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
	582	vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
[5133]	583
	584	! Also need ice/ocean stress on T points so that taum can be updated
	585	! This interpolation is already done in CICE so best to use those values
	586	CALL cice2nemo(strocnxT,ztmp1,'T',-1.)
	587	CALL cice2nemo(strocnyT,ztmp2,'T',-1.)
	588
	589	! Update taum with modulus of ice-ocean stress
	590	! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
[5836]	591	taum(:,:)=(1.0-fr_i(:,:))taum(:,:)+fr_i(:,:)SQRT(ztmp1ztmp1 + ztmp2ztmp2)
[2874]	592
	593	! Freshwater fluxes
	594
[4990]	595	IF (ksbc == jp_flx) THEN
[2874]	596	! Note that emp from the forcing files is evap(1-aice)-(tprecip-aicesprecip)
	597	! What we want here is evap(1-aice)-tprecip(1-aice) hence manipulation below
	598	! Not ideal since aice won't be the same as in the atmosphere.
	599	! Better to use evap and tprecip? (but for now don't read in evap in this case)
[3193]	600	emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
[7646]	601	ELSE IF (ksbc == jp_blk) THEN
[3193]	602	emp(:,:) = (1.0-fr_i(:,:))*emp(:,:)
[5407]	603	ELSE IF (ksbc == jp_purecpl) THEN
[3625]	604	! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
	605	! This is currently as required with the coupling fields from the UM atmosphere
[3193]	606	emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
	607	ENDIF
[2874]	608
[4990]	609	#if defined key_cice4
[3625]	610	CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
	611	CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
[4990]	612	#else
	613	CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
	614	CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
	615	#endif
[2874]	616
[3625]	617	! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
	618	! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
	619	! which has to be compensated for by the ocean salinity potentially going negative
	620	! This check breaks conservation but seems reasonable until we have prognostic ice salinity
	621	! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
	622	WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
	623	sfx(:,:)=ztmp2(:,:)*1000.0
	624	emp(:,:)=emp(:,:)-ztmp1(:,:)
[4990]	625	fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
	626
[9098]	627	CALL lbc_lnk_multi( emp , 'T', 1., sfx , 'T', 1. )
[2874]	628
	629	! Solar penetrative radiation and non solar surface heat flux
	630
	631	! Scale qsr and qns according to ice fraction (bulk formulae only)
	632
[7646]	633	IF (ksbc == jp_blk) THEN
[3193]	634	qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
	635	qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
	636	ENDIF
[2874]	637	! Take into account snow melting except for fully coupled when already in qns_tot
[5407]	638	IF (ksbc == jp_purecpl) THEN
[3193]	639	qsr(:,:)= qsr_tot(:,:)
	640	qns(:,:)= qns_tot(:,:)
	641	ELSE
	642	qns(:,:)= qns(:,:)-sprecip(:,:)Lfresh(1.0-fr_i(:,:))
	643	ENDIF
[2874]	644
	645	! Now add in ice / snow related terms
	646	! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
[4990]	647	#if defined key_cice4
[3625]	648	CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
[4990]	649	#else
	650	CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
	651	#endif
[3625]	652	qsr(:,:)=qsr(:,:)+ztmp1(:,:)
[3193]	653	CALL lbc_lnk( qsr , 'T', 1. )
[2874]	654
[3193]	655	DO jj=1,jpj
	656	DO ji=1,jpi
[2874]	657	nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
[3193]	658	ENDDO
	659	ENDDO
[2874]	660
[4990]	661	#if defined key_cice4
[3625]	662	CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
[4990]	663	#else
	664	CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
	665	#endif
[3625]	666	qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
[2874]	667
[3193]	668	CALL lbc_lnk( qns , 'T', 1. )
[2874]	669
	670	! Prepare for the following CICE time-step
	671
[3193]	672	CALL cice2nemo(aice,fr_i,'T', 1. )
[5407]	673	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
[3625]	674	DO jl=1,ncat
	675	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	676	ENDDO
	677	ENDIF
[2874]	678
	679	! T point to U point
	680	! T point to V point
[3193]	681	DO jj=1,jpjm1
	682	DO ji=1,jpim1
	683	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	684	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	685	ENDDO
	686	ENDDO
[2874]	687
[9098]	688	CALL lbc_lnk_multi( fr_iu , 'U', 1., fr_iv , 'V', 1. )
[2874]	689
[9019]	690	! set the snow+ice mass
	691	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	692	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	693	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	694	snwice_mass_b(:,:) = snwice_mass(:,:)
	695	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
[3193]	696	!
[2874]	697	END SUBROUTINE cice_sbc_out
	698
[3152]	699
[2874]	700	SUBROUTINE cice_sbc_hadgam( kt )
	701	!!---------------------------------------------------------------------
	702	!! * ROUTINE cice_sbc_hadgam *
[3040]	703	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
[2874]	704	!!
	705	!!
[9124]	706	!!---------------------------------------------------------------------
[2874]	707	INTEGER, INTENT( in ) :: kt ! ocean time step
[9124]	708	!!
[3625]	709	INTEGER :: jl ! dummy loop index
[3193]	710	INTEGER :: ierror
[9124]	711	!!---------------------------------------------------------------------
[9500]	712	#if defined key_asminc
	713	! Import fresh water and salt flux due to seaice da
	714	CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
	715	CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
	716	#endif
	717
[9498]	718	!
[2874]	719	! ! =========================== !
	720	! ! Prepare Coupling fields !
	721	! ! =========================== !
[9019]	722	!
	723	! x and y comp of ice velocity
	724	!
[3193]	725	CALL cice2nemo(uvel,u_ice,'F', -1. )
	726	CALL cice2nemo(vvel,v_ice,'F', -1. )
	727	!
[9019]	728	! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
	729	!
	730	! Snow and ice thicknesses (CO_2 and CO_3)
	731	!
	732	DO jl = 1, ncat
	733	CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
	734	CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
	735	END DO
[9500]	736
	737	#if ! defined key_cice4
	738	! Meltpond fraction and depth
	739	DO jl = 1,ncat
[9675]	740	CALL cice2nemo(apeffn(:,:,jl,:),a_ip(:,:,jl),'T', 1. )
	741	CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),v_ip(:,:,jl),'T', 1. )
[9500]	742	ENDDO
[9695]	743
	744	v_ip(:,:,:) = v_ip(:,:,:) * a_i(:,:,:)
[9500]	745	#endif
	746
	747
	748	! If using multilayers thermodynamics in CICE then get top layer temperature
	749	! and effective conductivity
	750	!! When using NEMO with CICE, this change requires use of
	751	!! one of the following two CICE branches:
	752	!! - at CICE5.0, hadax/r1015_GSI8_with_GSI7
	753	!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
	754	IF (heat_capacity) THEN
	755	DO jl = 1,ncat
[9675]	756	CALL cice2nemo(Tn_top(:,:,jl,:),t1_ice(:,:,jl),'T', 1. )
	757	CALL cice2nemo(keffn_top(:,:,jl,:),cnd_ice(:,:,jl),'T', 1. )
[9500]	758	ENDDO
	759	! Convert surface temperature to Kelvin
[9675]	760	t1_ice(:,:,:)=t1_ice(:,:,:)+rt0
[9500]	761	ELSE
[9675]	762	t1_ice(:,:,:) = 0.0
	763	cnd_ice(:,:,:) = 0.0
[9500]	764	ENDIF
	765
[9019]	766	!
[2874]	767	END SUBROUTINE cice_sbc_hadgam
	768
	769
	770	SUBROUTINE cice_sbc_final
	771	!!---------------------------------------------------------------------
	772	!! * ROUTINE cice_sbc_final *
	773	!! ** Purpose: Finalize CICE
	774	!!---------------------------------------------------------------------
[9124]	775	!
[2874]	776	IF(lwp) WRITE(numout,*)'cice_sbc_final'
[9124]	777	!
[3193]	778	CALL CICE_Finalize
[9124]	779	!
[2874]	780	END SUBROUTINE cice_sbc_final
	781
[9124]	782
[2874]	783	SUBROUTINE cice_sbc_force (kt)
	784	!!---------------------------------------------------------------------
	785	!! * ROUTINE cice_sbc_force *
	786	!! ** Purpose : Provide CICE forcing from files
	787	!!
	788	!!---------------------------------------------------------------------
	789	!! ** Method : READ monthly flux file in NetCDF files
	790	!!
	791	!! snowfall
	792	!! rainfall
	793	!! sublimation rate
	794	!! topmelt (category)
	795	!! botmelt (category)
	796	!!
	797	!! History :
	798	!!----------------------------------------------------------------------
	799	USE iom
[9124]	800	!!
	801	INTEGER, INTENT( in ) :: kt ! ocean time step
	802	!!
[2874]	803	INTEGER :: ierror ! return error code
	804	INTEGER :: ifpr ! dummy loop index
	805	!!
	806	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
	807	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
	808	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
	809	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
	810	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
	811	!!
	812	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
	813	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
	814	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
[4230]	815	INTEGER :: ios
[2874]	816	!!---------------------------------------------------------------------
	817
	818	! ! ====================== !
	819	IF( kt == nit000 ) THEN ! First call kt=nit000 !
	820	! ! ====================== !
[4990]	821	! namsbc_cice is not yet in the reference namelist
	822	! set file information (default values)
	823	cn_dir = './' ! directory in which the model is executed
	824
	825	! (NB: frequency positive => hours, negative => months)
	826	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask
	827	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file
	828	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	829	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	830	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' )
	831	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	832	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	833	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	834	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	835	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	836	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	837	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	838	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	839	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	840	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	841
[4230]	842	REWIND( numnam_ref ) ! Namelist namsbc_cice in reference namelist :
	843	READ ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
[9168]	844	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
[2874]	845
[4230]	846	REWIND( numnam_cfg ) ! Namelist namsbc_cice in configuration namelist : Parameters of the run
	847	READ ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
[9168]	848	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
[4624]	849	IF(lwm) WRITE ( numond, namsbc_cice )
[2874]	850
	851	! store namelist information in an array
	852	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
	853	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
	854	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
	855	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
	856	slf_i(jp_bot5) = sn_bot5
	857
	858	! set sf structure
	859	ALLOCATE( sf(jpfld), STAT=ierror )
	860	IF( ierror > 0 ) THEN
	861	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
	862	ENDIF
	863
	864	DO ifpr= 1, jpfld
	865	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
	866	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
	867	END DO
	868
	869	! fill sf with slf_i and control print
	870	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
	871	!
	872	ENDIF
	873
	874	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
	875	! ! input fields at the current time-step
	876
	877	! set the fluxes from read fields
	878	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
	879	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
[3040]	880	! May be better to do this conversion somewhere else
[9675]	881	zevap_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
	882	qml_ice(:,:,1) = sf(jp_top1)%fnow(:,:,1)
	883	qml_ice(:,:,2) = sf(jp_top2)%fnow(:,:,1)
	884	qml_ice(:,:,3) = sf(jp_top3)%fnow(:,:,1)
	885	qml_ice(:,:,4) = sf(jp_top4)%fnow(:,:,1)
	886	qml_ice(:,:,5) = sf(jp_top5)%fnow(:,:,1)
	887	qcn_ice(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
	888	qcn_ice(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
	889	qcn_ice(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
	890	qcn_ice(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
	891	qcn_ice(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
[2874]	892
	893	! control print (if less than 100 time-step asked)
	894	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
	895	WRITE(numout,*)
	896	WRITE(numout,*) ' read forcing fluxes for CICE OK'
	897	CALL FLUSH(numout)
	898	ENDIF
	899
	900	END SUBROUTINE cice_sbc_force
	901
	902	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
	903	!!---------------------------------------------------------------------
	904	!! * ROUTINE nemo2cice *
	905	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
	906	#if defined key_nemocice_decomp
	907	!!
	908	!! NEMO and CICE PE sub domains are identical, hence
	909	!! there is no need to gather or scatter data from
	910	!! one PE configuration to another.
	911	#else
	912	!! Automatically gather/scatter between
	913	!! different processors and blocks
	914	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
	915	!! B. Gather pn into global array (png)
	916	!! C. Map png into CICE global array (pcg)
	917	!! D. Scatter pcg to CICE blocks (pc) + update haloes
	918	#endif
	919	!!---------------------------------------------------------------------
[3193]	920	CHARACTER(len=1), INTENT( in ) :: &
	921	cd_type ! nature of pn grid-point
	922	! ! = T or F gridpoints
	923	REAL(wp), INTENT( in ) :: &
	924	psgn ! control of the sign change
	925	! ! =-1 , the sign is modified following the type of b.c. used
	926	! ! = 1 , no sign change
	927	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	928	#if !defined key_nemocice_decomp
[3625]	929	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
[3193]	930	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	931	#endif
[3193]	932	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
	933	INTEGER (int_kind) :: &
	934	field_type, &! id for type of field (scalar, vector, angle)
	935	grid_loc ! id for location on horizontal grid
[2874]	936	! (center, NEcorner, Nface, Eface)
	937
[3193]	938	INTEGER :: ji, jj, jn ! dummy loop indices
[5836]	939	!!---------------------------------------------------------------------
[2874]	940
[3193]	941	! A. Ensure all haloes are filled in NEMO field (pn)
[2874]	942
[3193]	943	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	944
	945	#if defined key_nemocice_decomp
	946
[3193]	947	! Copy local domain data from NEMO to CICE field
	948	pc(:,:,1)=0.0
[9631]	949	DO jj=2,jpj-1
	950	DO ji=2,jpi-1
[3625]	951	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
[3193]	952	ENDDO
	953	ENDDO
[2874]	954
	955	#else
	956
[3193]	957	! B. Gather pn into global array (png)
[2874]	958
[3193]	959	IF ( jpnij > 1) THEN
	960	CALL mppsync
	961	CALL mppgather (pn,0,png)
	962	CALL mppsync
	963	ELSE
	964	png(:,:,1)=pn(:,:)
	965	ENDIF
[2874]	966
[3193]	967	! C. Map png into CICE global array (pcg)
[2874]	968
	969	! Need to make sure this is robust to changes in NEMO halo rows....
	970	! (may be OK but not 100% sure)
	971
[3193]	972	IF (nproc==0) THEN
[2874]	973	! pcg(:,:)=0.0
[3193]	974	DO jn=1,jpnij
[3625]	975	DO jj=nldjt(jn),nlejt(jn)
	976	DO ji=nldit(jn),nleit(jn)
	977	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
[3193]	978	ENDDO
	979	ENDDO
	980	ENDDO
[3625]	981	DO jj=1,ny_global
	982	DO ji=1,nx_global
	983	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
	984	ENDDO
	985	ENDDO
[3193]	986	ENDIF
[2874]	987
	988	#endif
	989
[3193]	990	SELECT CASE ( cd_type )
	991	CASE ( 'T' )
	992	grid_loc=field_loc_center
	993	CASE ( 'F' )
	994	grid_loc=field_loc_NEcorner
	995	END SELECT
[2874]	996
[3193]	997	SELECT CASE ( NINT(psgn) )
	998	CASE ( -1 )
	999	field_type=field_type_vector
	1000	CASE ( 1 )
	1001	field_type=field_type_scalar
	1002	END SELECT
[2874]	1003
	1004	#if defined key_nemocice_decomp
[3193]	1005	! Ensure CICE halos are up to date
	1006	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1007	#else
[3193]	1008	! D. Scatter pcg to CICE blocks (pc) + update halos
	1009	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
[2874]	1010	#endif
	1011
	1012	END SUBROUTINE nemo2cice
	1013
	1014	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
	1015	!!---------------------------------------------------------------------
	1016	!! * ROUTINE cice2nemo *
	1017	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
	1018	#if defined key_nemocice_decomp
	1019	!!
	1020	!! NEMO and CICE PE sub domains are identical, hence
	1021	!! there is no need to gather or scatter data from
	1022	!! one PE configuration to another.
	1023	#else
	1024	!! Automatically deal with scatter/gather between
	1025	!! different processors and blocks
	1026	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
	1027	!! B. Map pcg into NEMO global array (png)
	1028	!! C. Scatter png into NEMO field (pn) for each processor
	1029	!! D. Ensure all haloes are filled in pn
	1030	#endif
	1031	!!---------------------------------------------------------------------
	1032
[3193]	1033	CHARACTER(len=1), INTENT( in ) :: &
	1034	cd_type ! nature of pn grid-point
	1035	! ! = T or F gridpoints
	1036	REAL(wp), INTENT( in ) :: &
	1037	psgn ! control of the sign change
	1038	! ! =-1 , the sign is modified following the type of b.c. used
	1039	! ! = 1 , no sign change
	1040	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	1041
	1042	#if defined key_nemocice_decomp
[3193]	1043	INTEGER (int_kind) :: &
	1044	field_type, & ! id for type of field (scalar, vector, angle)
	1045	grid_loc ! id for location on horizontal grid
	1046	! (center, NEcorner, Nface, Eface)
[2874]	1047	#else
[3193]	1048	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	1049	#endif
	1050
[3193]	1051	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
[2874]	1052
[3193]	1053	INTEGER :: ji, jj, jn ! dummy loop indices
[2874]	1054
	1055
	1056	#if defined key_nemocice_decomp
	1057
[3193]	1058	SELECT CASE ( cd_type )
	1059	CASE ( 'T' )
	1060	grid_loc=field_loc_center
	1061	CASE ( 'F' )
	1062	grid_loc=field_loc_NEcorner
	1063	END SELECT
[2874]	1064
[3193]	1065	SELECT CASE ( NINT(psgn) )
	1066	CASE ( -1 )
	1067	field_type=field_type_vector
	1068	CASE ( 1 )
	1069	field_type=field_type_scalar
	1070	END SELECT
[2874]	1071
[3193]	1072	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1073
	1074
[3193]	1075	pn(:,:)=0.0
	1076	DO jj=1,jpjm1
	1077	DO ji=1,jpim1
[3625]	1078	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
[3193]	1079	ENDDO
	1080	ENDDO
[2874]	1081
	1082	#else
	1083
[3193]	1084	! A. Gather CICE blocks (pc) into global array (pcg)
[2874]	1085
[3193]	1086	CALL gather_global(pcg, pc, 0, distrb_info)
[2874]	1087
	1088	! B. Map pcg into NEMO global array (png)
	1089
	1090	! Need to make sure this is robust to changes in NEMO halo rows....
	1091	! (may be OK but not spent much time thinking about it)
[3625]	1092	! Note that non-existent pcg elements may be used below, but
	1093	! the lbclnk call on pn will replace these with sensible values
[2874]	1094
[3193]	1095	IF (nproc==0) THEN
	1096	png(:,:,:)=0.0
	1097	DO jn=1,jpnij
[3625]	1098	DO jj=nldjt(jn),nlejt(jn)
	1099	DO ji=nldit(jn),nleit(jn)
	1100	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
[3193]	1101	ENDDO
	1102	ENDDO
	1103	ENDDO
	1104	ENDIF
[2874]	1105
[3193]	1106	! C. Scatter png into NEMO field (pn) for each processor
[2874]	1107
[3193]	1108	IF ( jpnij > 1) THEN
	1109	CALL mppsync
	1110	CALL mppscatter (png,0,pn)
	1111	CALL mppsync
	1112	ELSE
	1113	pn(:,:)=png(:,:,1)
	1114	ENDIF
[2874]	1115
	1116	#endif
	1117
[3193]	1118	! D. Ensure all haloes are filled in pn
[2874]	1119
[3193]	1120	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	1121
	1122	END SUBROUTINE cice2nemo
	1123
	1124	#else
	1125	!!----------------------------------------------------------------------
	1126	!! Default option Dummy module NO CICE sea-ice model
	1127	!!----------------------------------------------------------------------
	1128	CONTAINS
	1129
[4990]	1130	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
[2874]	1131	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
	1132	END SUBROUTINE sbc_ice_cice
	1133
[4990]	1134	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
[2874]	1135	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
	1136	END SUBROUTINE cice_sbc_init
	1137
	1138	SUBROUTINE cice_sbc_final ! Dummy routine
	1139	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
	1140	END SUBROUTINE cice_sbc_final
	1141
	1142	#endif
	1143
	1144	!!======================================================================
	1145	END MODULE sbcice_cice

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

Download in other formats: