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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_w_keyasminc/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 7749

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

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

Download in other formats: