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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8037

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

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

Download in other formats: