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 NEMO/branches/2020/ticket2487/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/ticket2487/src/OCE/SBC/sbcice_cice.F90 @ 13179

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

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

Download in other formats: