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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_fix_zemp_ice_10681/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11029

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

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

Download in other formats: