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_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5796

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

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

Download in other formats: