Context Navigation

sbcice_cice.F90 @ 8933

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8933

Download in other formats: