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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11935

Last change on this file since 11935 was 11935, checked in by dcarneir, 4 years ago
Changing SBC scripts (sbc_oce and sbcice_cice) to include sea ice thickness
File size: 45.7 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
[3625]	17	USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic
[2874]	18	USE in_out_manager ! I/O manager
[4990]	19	USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit
[2874]	20	USE lib_mpp ! distributed memory computing library
	21	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[3186]	22	USE wrk_nemo ! work arrays
[3193]	23	USE timing ! Timing
[2874]	24	USE daymod ! calendar
	25	USE fldread ! read input fields
	26	USE sbc_oce ! Surface boundary condition: ocean fields
	27	USE sbc_ice ! Surface boundary condition: ice fields
	28	USE sbcblk_core ! Surface boundary condition: CORE bulk
	29	USE sbccpl
	30
	31	USE ice_kinds_mod
	32	USE ice_blocks
	33	USE ice_domain
	34	USE ice_domain_size
	35	USE ice_boundary
	36	USE ice_constants
	37	USE ice_gather_scatter
	38	USE ice_calendar, only: dt
[3625]	39	USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
[4990]	40	# if defined key_cice4
[2874]	41	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[5133]	42	strocnxT,strocnyT, &
[3189]	43	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, &
	44	fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, &
[2874]	45	flatn_f,fsurfn_f,fcondtopn_f, &
	46	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
	47	swvdr,swvdf,swidr,swidf
[4990]	48	USE ice_therm_vertical, only: calc_Tsfc
	49	#else
	50	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[5133]	51	strocnxT,strocnyT, &
[4990]	52	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai, &
	53	fresh_ai,fhocn_ai,fswthru_ai,frzmlt, &
	54	flatn_f,fsurfn_f,fcondtopn_f, &
	55	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
	56	swvdr,swvdf,swidr,swidf
	57	USE ice_therm_shared, only: calc_Tsfc
	58	#endif
[2874]	59	USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
[3176]	60	USE ice_atmo, only: calc_strair
[2874]	61
	62	USE CICE_InitMod
	63	USE CICE_RunMod
	64	USE CICE_FinalMod
	65
	66	IMPLICIT NONE
	67	PRIVATE
	68
	69	!! * Routine accessibility
	70	PUBLIC cice_sbc_init ! routine called by sbc_init
	71	PUBLIC cice_sbc_final ! routine called by sbc_final
	72	PUBLIC sbc_ice_cice ! routine called by sbc
	73
[4627]	74	INTEGER :: ji_off
	75	INTEGER :: jj_off
[3625]	76
[2874]	77	INTEGER , PARAMETER :: jpfld = 13 ! maximum number of files to read
	78	INTEGER , PARAMETER :: jp_snow = 1 ! index of snow file
	79	INTEGER , PARAMETER :: jp_rain = 2 ! index of rain file
	80	INTEGER , PARAMETER :: jp_sblm = 3 ! index of sublimation file
	81	INTEGER , PARAMETER :: jp_top1 = 4 ! index of category 1 topmelt file
	82	INTEGER , PARAMETER :: jp_top2 = 5 ! index of category 2 topmelt file
	83	INTEGER , PARAMETER :: jp_top3 = 6 ! index of category 3 topmelt file
	84	INTEGER , PARAMETER :: jp_top4 = 7 ! index of category 4 topmelt file
	85	INTEGER , PARAMETER :: jp_top5 = 8 ! index of category 5 topmelt file
	86	INTEGER , PARAMETER :: jp_bot1 = 9 ! index of category 1 botmelt file
	87	INTEGER , PARAMETER :: jp_bot2 = 10 ! index of category 2 botmelt file
	88	INTEGER , PARAMETER :: jp_bot3 = 11 ! index of category 3 botmelt file
	89	INTEGER , PARAMETER :: jp_bot4 = 12 ! index of category 4 botmelt file
	90	INTEGER , PARAMETER :: jp_bot5 = 13 ! index of category 5 botmelt file
	91	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
	92
	93	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE :: png ! local array used in sbc_cice_ice
	94
	95	!! * Substitutions
	96	# include "domzgr_substitute.h90"
	97
[5215]	98	!! $Id$
[2874]	99	CONTAINS
	100
	101	INTEGER FUNCTION sbc_ice_cice_alloc()
	102	!!----------------------------------------------------------------------
	103	!! * FUNCTION sbc_ice_cice_alloc *
	104	!!----------------------------------------------------------------------
	105	ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
	106	IF( lk_mpp ) CALL mpp_sum ( sbc_ice_cice_alloc )
	107	IF( sbc_ice_cice_alloc > 0 ) CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
	108	END FUNCTION sbc_ice_cice_alloc
	109
[4990]	110	SUBROUTINE sbc_ice_cice( kt, ksbc )
[2874]	111	!!---------------------------------------------------------------------
	112	!! * ROUTINE sbc_ice_cice *
	113	!!
	114	!! ** Purpose : update the ocean surface boundary condition via the
	115	!! CICE Sea Ice Model time stepping
	116	!!
[3040]	117	!! ** Method : - Get any extra forcing fields for CICE
	118	!! - Prepare forcing fields
[2874]	119	!! - CICE model time stepping
	120	!! - call the routine that computes mass and
	121	!! heat fluxes at the ice/ocean interface
	122	!!
	123	!! ** Action : - time evolution of the CICE sea-ice model
	124	!! - update all sbc variables below sea-ice:
[3625]	125	!! utau, vtau, qns , qsr, emp , sfx
[2874]	126	!!---------------------------------------------------------------------
	127	INTEGER, INTENT(in) :: kt ! ocean time step
[4990]	128	INTEGER, INTENT(in) :: ksbc ! surface forcing type
[2874]	129	!!----------------------------------------------------------------------
[3193]	130	!
	131	IF( nn_timing == 1 ) CALL timing_start('sbc_ice_cice')
	132	!
[2874]	133	! !----------------------!
	134	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Ice time-step only !
	135	! !----------------------!
	136
	137	! Make sure any fluxes required for CICE are set
[4990]	138	IF ( ksbc == jp_flx ) THEN
[2874]	139	CALL cice_sbc_force(kt)
[5407]	140	ELSE IF ( ksbc == jp_purecpl ) THEN
[2874]	141	CALL sbc_cpl_ice_flx( 1.0-fr_i )
	142	ENDIF
	143
[4990]	144	CALL cice_sbc_in ( kt, ksbc )
[2874]	145	CALL CICE_Run
[4990]	146	CALL cice_sbc_out ( kt, ksbc )
[2874]	147
[5407]	148	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(kt+1)
[2874]	149
	150	ENDIF ! End sea-ice time step only
[3193]	151	!
	152	IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_cice')
[2874]	153
	154	END SUBROUTINE sbc_ice_cice
	155
[4990]	156	SUBROUTINE cice_sbc_init (ksbc)
[2874]	157	!!---------------------------------------------------------------------
	158	!! * ROUTINE cice_sbc_init *
[3040]	159	!! ** Purpose: Initialise ice related fields for NEMO and coupling
[2874]	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
[4990]	193	ELSEIF (ksbc == jp_core) 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
[11935]	232	! Snow and ice thickness
	233	! vice, vsno is volume per unit area of grid cell = thickness
	234	CALL cice2nemo(vice,thick_i,'T', 1. )
	235	CALL cice2nemo(vsno,thick_s,'T', 1. )
	236
	237	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
	238	DO jl = 1,ncat
	239	CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1.)
	240	CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1.)
	241	ENDDO
	242	ENDIF
	243
	244	! T point to U point
	245	! T point to V point
	246
	247	! Sea ice thickness
	248	thick_iu(:,:)=0.0
	249	thick_iv(:,:)=0.0
	250	DO jj=1,jpjm1
	251	DO ji=1,jpim1
	252	thick_iu(ji,jj)=0.5(thick_i(ji,jj)+thick_i(ji+1,jj))umask(ji,jj,1)
	253	thick_iv(ji,jj)=0.5(thick_i(ji,jj)+thick_i(ji,jj+1))vmask(ji,jj,1)
	254	ENDDO
	255	ENDDO
	256
	257	CALL lbc_lnk ( thick_iu , 'U', 1. )
	258	CALL lbc_lnk ( thick_iv , 'V', 1. )
	259
	260	! Snow depth
	261	thick_su(:,:)=0.0
	262	thick_sv(:,:)=0.0
	263	DO jj=1,jpjm1
	264	DO ji=1,jpim1
	265	thick_su(ji,jj)=0.5(thick_s(ji,jj)+thick_s(ji+1,jj))umask(ji,jj,1)
	266	thick_sv(ji,jj)=0.5(thick_s(ji,jj)+thick_s(ji,jj+1))vmask(ji,jj,1)
	267	ENDDO
	268	ENDDO
	269
	270	CALL lbc_lnk ( thick_su , 'U', 1. )
	271	CALL lbc_lnk ( thick_sv , 'V', 1. )
	272
[3625]	273	! ! embedded sea ice
	274	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
	275	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	276	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	277	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	278	snwice_mass_b(:,:) = snwice_mass(:,:)
	279	ELSE
	280	snwice_mass (:,:) = 0.0_wp ! no mass exchanges
	281	snwice_mass_b(:,:) = 0.0_wp ! no mass exchanges
	282	ENDIF
[4990]	283	IF( .NOT. ln_rstart ) THEN
	284	IF( nn_ice_embd == 2 ) THEN ! full embedment (case 2) deplete the initial ssh below sea-ice area
	285	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	286	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
	287	#if defined key_vvl
	288	! key_vvl necessary? clem: yes for compilation purpose
	289	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	290	fse3t_n(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshn(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
	291	fse3t_b(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshb(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
	292	ENDDO
	293	fse3t_a(:,:,:) = fse3t_b(:,:,:)
	294	! Reconstruction of all vertical scale factors at now and before time
	295	! steps
	296	! =============================================================================
	297	! Horizontal scale factor interpolations
	298	! --------------------------------------
	299	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
	300	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
	301	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
	302	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
	303	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
	304	! Vertical scale factor interpolations
	305	! ------------------------------------
	306	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
	307	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
	308	CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
	309	CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
	310	CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
	311	! t- and w- points depth
	312	! ----------------------
	313	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
	314	fsdepw_n(:,:,1) = 0.0_wp
	315	fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
	316	DO jk = 2, jpk
	317	fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
	318	fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
	319	fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:)
	320	END DO
	321	#endif
	322	ENDIF
[3625]	323	ENDIF
	324
	325	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
[3193]	326	!
	327	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_init')
	328	!
[2874]	329	END SUBROUTINE cice_sbc_init
	330
[3152]	331
[4990]	332	SUBROUTINE cice_sbc_in (kt, ksbc)
[2874]	333	!!---------------------------------------------------------------------
	334	!! * ROUTINE cice_sbc_in *
[3040]	335	!! ** Purpose: Set coupling fields and pass to CICE
[2874]	336	!!---------------------------------------------------------------------
[3152]	337	INTEGER, INTENT(in ) :: kt ! ocean time step
[4990]	338	INTEGER, INTENT(in ) :: ksbc ! surface forcing type
[2874]	339
[3625]	340	INTEGER :: ji, jj, jl ! dummy loop indices
	341	REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice
[3152]	342	REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
[3625]	343	REAL(wp) :: zintb, zintn ! dummy argument
[3152]	344	!!---------------------------------------------------------------------
[2874]	345
[3193]	346	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_in')
	347	!
[3625]	348	CALL wrk_alloc( jpi,jpj, ztmp, zpice )
[3152]	349	CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
[2874]	350
[3193]	351	IF( kt == nit000 ) THEN
[2874]	352	IF(lwp) WRITE(numout,*)'cice_sbc_in'
[3193]	353	ENDIF
[2874]	354
[3193]	355	ztmp(:,:)=0.0
[2874]	356
	357	! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
	358	! the first time-step)
	359
	360	! forced and coupled case
	361
[5407]	362	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
[2874]	363
[3193]	364	ztmpn(:,:,:)=0.0
[2874]	365
	366	! x comp of wind stress (CI_1)
	367	! U point to F point
[3193]	368	DO jj=1,jpjm1
	369	DO ji=1,jpi
	370	ztmp(ji,jj) = 0.5 * ( fr_iu(ji,jj) * utau(ji,jj) &
	371	+ fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
	372	ENDDO
	373	ENDDO
	374	CALL nemo2cice(ztmp,strax,'F', -1. )
[2874]	375
	376	! y comp of wind stress (CI_2)
	377	! V point to F point
[3193]	378	DO jj=1,jpj
	379	DO ji=1,jpim1
	380	ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
	381	+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
	382	ENDDO
	383	ENDDO
	384	CALL nemo2cice(ztmp,stray,'F', -1. )
[2874]	385
	386	! Surface downward latent heat flux (CI_5)
[4990]	387	IF (ksbc == jp_flx) THEN
[3625]	388	DO jl=1,ncat
	389	ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
[3193]	390	ENDDO
	391	ELSE
[2874]	392	! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
[3193]	393	qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub
[2874]	394	! End of temporary code
[3193]	395	DO jj=1,jpj
	396	DO ji=1,jpi
	397	IF (fr_i(ji,jj).eq.0.0) THEN
[3625]	398	DO jl=1,ncat
	399	ztmpn(ji,jj,jl)=0.0
[3193]	400	ENDDO
	401	! This will then be conserved in CICE
	402	ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
	403	ELSE
[3625]	404	DO jl=1,ncat
	405	ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
[3193]	406	ENDDO
	407	ENDIF
	408	ENDDO
	409	ENDDO
	410	ENDIF
[3625]	411	DO jl=1,ncat
	412	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
[2874]	413
	414	! GBM conductive flux through ice (CI_6)
	415	! Convert to GBM
[4990]	416	IF (ksbc == jp_flx) THEN
[3625]	417	ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
[3193]	418	ELSE
[3625]	419	ztmp(:,:) = botmelt(:,:,jl)
[3193]	420	ENDIF
[3625]	421	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
[2874]	422
	423	! GBM surface heat flux (CI_7)
	424	! Convert to GBM
[4990]	425	IF (ksbc == jp_flx) THEN
[3625]	426	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)
[3193]	427	ELSE
[3625]	428	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
[3193]	429	ENDIF
[3625]	430	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
[3193]	431	ENDDO
[2874]	432
[4990]	433	ELSE IF (ksbc == jp_core) THEN
[2874]	434
	435	! Pass CORE forcing fields to CICE (which will calculate heat fluxes etc itself)
	436	! x comp and y comp of atmosphere surface wind (CICE expects on T points)
[3193]	437	ztmp(:,:) = wndi_ice(:,:)
	438	CALL nemo2cice(ztmp,uatm,'T', -1. )
	439	ztmp(:,:) = wndj_ice(:,:)
	440	CALL nemo2cice(ztmp,vatm,'T', -1. )
	441	ztmp(:,:) = SQRT ( wndi_ice(:,:)2 + wndj_ice(:,:)2 )
	442	CALL nemo2cice(ztmp,wind,'T', 1. ) ! Wind speed (m/s)
	443	ztmp(:,:) = qsr_ice(:,:,1)
	444	CALL nemo2cice(ztmp,fsw,'T', 1. ) ! Incoming short-wave (W/m^2)
	445	ztmp(:,:) = qlw_ice(:,:,1)
	446	CALL nemo2cice(ztmp,flw,'T', 1. ) ! Incoming long-wave (W/m^2)
	447	ztmp(:,:) = tatm_ice(:,:)
	448	CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K)
	449	CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K)
[2874]	450	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
[3193]	451	ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
	452	! Constant (101000.) atm pressure assumed
	453	CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3)
	454	ztmp(:,:) = qatm_ice(:,:)
	455	CALL nemo2cice(ztmp,Qa,'T', 1. ) ! Specific humidity (kg/kg)
	456	ztmp(:,:)=10.0
	457	CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m)
[2874]	458
	459	! May want to check all values are physically realistic (as in CICE routine
	460	! prepare_forcing)?
	461
	462	! Divide shortwave into spectral bands (as in prepare_forcing)
[3193]	463	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
[2874]	464	CALL nemo2cice(ztmp,swvdr,'T', 1. )
[3193]	465	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
[2874]	466	CALL nemo2cice(ztmp,swvdf,'T', 1. )
[3193]	467	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
[2874]	468	CALL nemo2cice(ztmp,swidr,'T', 1. )
[3193]	469	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
[2874]	470	CALL nemo2cice(ztmp,swidf,'T', 1. )
	471
	472	ENDIF
	473
	474	! Snowfall
[4990]	475	! Ensure fsnow is positive (as in CICE routine prepare_forcing)
	476	IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
[3193]	477	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
	478	CALL nemo2cice(ztmp,fsnow,'T', 1. )
[2874]	479
	480	! Rainfall
[4990]	481	IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
[3193]	482	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
	483	CALL nemo2cice(ztmp,frain,'T', 1. )
[2874]	484
	485	! Freezing/melting potential
[3275]	486	! Calculated over NEMO leapfrog timestep (hence 2*dt)
[4292]	487	nfrzmlt(:,:)=rau0rcpfse3t_m(:,:)(Tocnfrz-sst_m(:,:))/(2.0dt)
[2874]	488
[3193]	489	ztmp(:,:) = nfrzmlt(:,:)
	490	CALL nemo2cice(ztmp,frzmlt,'T', 1. )
[2874]	491
	492	! SST and SSS
	493
[3193]	494	CALL nemo2cice(sst_m,sst,'T', 1. )
	495	CALL nemo2cice(sss_m,sss,'T', 1. )
[2874]	496
	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
[2874]	728	! Release work space
	729
[3625]	730	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
[3193]	731	!
	732	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_out')
	733	!
[2874]	734	END SUBROUTINE cice_sbc_out
	735
[3152]	736
[2874]	737	SUBROUTINE cice_sbc_hadgam( kt )
	738	!!---------------------------------------------------------------------
	739	!! * ROUTINE cice_sbc_hadgam *
[3040]	740	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
[2874]	741	!!
	742	!!
	743	INTEGER, INTENT( in ) :: kt ! ocean time step
	744	!!---------------------------------------------------------------------
	745
[3625]	746	INTEGER :: jl ! dummy loop index
[3193]	747	INTEGER :: ierror
[2874]	748
[3193]	749	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_hadgam')
	750	!
	751	IF( kt == nit000 ) THEN
[2874]	752	IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
	753	IF( sbc_cpl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
[3193]	754	ENDIF
[2874]	755
	756	! ! =========================== !
	757	! ! Prepare Coupling fields !
	758	! ! =========================== !
	759
	760	! x and y comp of ice velocity
	761
[3193]	762	CALL cice2nemo(uvel,u_ice,'F', -1. )
	763	CALL cice2nemo(vvel,v_ice,'F', -1. )
[2874]	764
	765	! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
	766
	767	! Snow and ice thicknesses (CO_2 and CO_3)
	768
[3625]	769	DO jl = 1,ncat
	770	CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. )
	771	CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
[3193]	772	ENDDO
	773	!
	774	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_hadgam')
	775	!
[2874]	776	END SUBROUTINE cice_sbc_hadgam
	777
	778
	779	SUBROUTINE cice_sbc_final
	780	!!---------------------------------------------------------------------
	781	!! * ROUTINE cice_sbc_final *
	782	!! ** Purpose: Finalize CICE
	783	!!---------------------------------------------------------------------
	784
	785	IF(lwp) WRITE(numout,*)'cice_sbc_final'
	786
[3193]	787	CALL CICE_Finalize
[2874]	788
	789	END SUBROUTINE cice_sbc_final
	790
	791	SUBROUTINE cice_sbc_force (kt)
	792	!!---------------------------------------------------------------------
	793	!! * ROUTINE cice_sbc_force *
	794	!! ** Purpose : Provide CICE forcing from files
	795	!!
	796	!!---------------------------------------------------------------------
	797	!! ** Method : READ monthly flux file in NetCDF files
	798	!!
	799	!! snowfall
	800	!! rainfall
	801	!! sublimation rate
	802	!! topmelt (category)
	803	!! botmelt (category)
	804	!!
	805	!! History :
	806	!!----------------------------------------------------------------------
	807	!! * Modules used
	808	USE iom
	809
	810	!! * arguments
	811	INTEGER, INTENT( in ) :: kt ! ocean time step
	812
	813	INTEGER :: ierror ! return error code
	814	INTEGER :: ifpr ! dummy loop index
	815	!!
	816	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
	817	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
	818	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
	819	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
	820	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
	821
	822	!!
	823	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
	824	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
	825	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
[4230]	826	INTEGER :: ios
[2874]	827	!!---------------------------------------------------------------------
	828
	829	! ! ====================== !
	830	IF( kt == nit000 ) THEN ! First call kt=nit000 !
	831	! ! ====================== !
[4990]	832	! namsbc_cice is not yet in the reference namelist
	833	! set file information (default values)
	834	cn_dir = './' ! directory in which the model is executed
	835
	836	! (NB: frequency positive => hours, negative => months)
	837	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask
	838	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file
	839	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	840	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	841	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' )
	842	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	843	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	844	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	845	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	846	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	847	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	848	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	849	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	850	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	851	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	852
[4230]	853	REWIND( numnam_ref ) ! Namelist namsbc_cice in reference namelist :
	854	READ ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
	855	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
[2874]	856
[4230]	857	REWIND( numnam_cfg ) ! Namelist namsbc_cice in configuration namelist : Parameters of the run
	858	READ ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
	859	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
[4624]	860	IF(lwm) WRITE ( numond, namsbc_cice )
[2874]	861
	862	! store namelist information in an array
	863	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
	864	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
	865	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
	866	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
	867	slf_i(jp_bot5) = sn_bot5
	868
	869	! set sf structure
	870	ALLOCATE( sf(jpfld), STAT=ierror )
	871	IF( ierror > 0 ) THEN
	872	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
	873	ENDIF
	874
	875	DO ifpr= 1, jpfld
	876	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
	877	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
	878	END DO
	879
	880	! fill sf with slf_i and control print
	881	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
	882	!
	883	ENDIF
	884
	885	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
	886	! ! input fields at the current time-step
	887
	888	! set the fluxes from read fields
	889	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
	890	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
[3040]	891	! May be better to do this conversion somewhere else
[2874]	892	qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
	893	topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
	894	topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
	895	topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
	896	topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
	897	topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
	898	botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
	899	botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
	900	botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
	901	botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
	902	botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
	903
	904	! control print (if less than 100 time-step asked)
	905	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
	906	WRITE(numout,*)
	907	WRITE(numout,*) ' read forcing fluxes for CICE OK'
	908	CALL FLUSH(numout)
	909	ENDIF
	910
	911	END SUBROUTINE cice_sbc_force
	912
	913	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
	914	!!---------------------------------------------------------------------
	915	!! * ROUTINE nemo2cice *
	916	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
	917	#if defined key_nemocice_decomp
	918	!!
	919	!! NEMO and CICE PE sub domains are identical, hence
	920	!! there is no need to gather or scatter data from
	921	!! one PE configuration to another.
	922	#else
	923	!! Automatically gather/scatter between
	924	!! different processors and blocks
	925	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
	926	!! B. Gather pn into global array (png)
	927	!! C. Map png into CICE global array (pcg)
	928	!! D. Scatter pcg to CICE blocks (pc) + update haloes
	929	#endif
	930	!!---------------------------------------------------------------------
	931
[3193]	932	CHARACTER(len=1), INTENT( in ) :: &
	933	cd_type ! nature of pn grid-point
	934	! ! = T or F gridpoints
	935	REAL(wp), INTENT( in ) :: &
	936	psgn ! control of the sign change
	937	! ! =-1 , the sign is modified following the type of b.c. used
	938	! ! = 1 , no sign change
	939	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	940	#if !defined key_nemocice_decomp
[3625]	941	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
[3193]	942	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	943	#endif
[3193]	944	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
	945	INTEGER (int_kind) :: &
	946	field_type, &! id for type of field (scalar, vector, angle)
	947	grid_loc ! id for location on horizontal grid
[2874]	948	! (center, NEcorner, Nface, Eface)
	949
[3193]	950	INTEGER :: ji, jj, jn ! dummy loop indices
[2874]	951
[3193]	952	! A. Ensure all haloes are filled in NEMO field (pn)
[2874]	953
[3193]	954	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	955
	956	#if defined key_nemocice_decomp
	957
[3193]	958	! Copy local domain data from NEMO to CICE field
	959	pc(:,:,1)=0.0
[3625]	960	DO jj=2,ny_block-1
	961	DO ji=2,nx_block-1
	962	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
[3193]	963	ENDDO
	964	ENDDO
[2874]	965
	966	#else
	967
[3193]	968	! B. Gather pn into global array (png)
[2874]	969
[3193]	970	IF ( jpnij > 1) THEN
	971	CALL mppsync
	972	CALL mppgather (pn,0,png)
	973	CALL mppsync
	974	ELSE
	975	png(:,:,1)=pn(:,:)
	976	ENDIF
[2874]	977
[3193]	978	! C. Map png into CICE global array (pcg)
[2874]	979
	980	! Need to make sure this is robust to changes in NEMO halo rows....
	981	! (may be OK but not 100% sure)
	982
[3193]	983	IF (nproc==0) THEN
[2874]	984	! pcg(:,:)=0.0
[3193]	985	DO jn=1,jpnij
[3625]	986	DO jj=nldjt(jn),nlejt(jn)
	987	DO ji=nldit(jn),nleit(jn)
	988	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
[3193]	989	ENDDO
	990	ENDDO
	991	ENDDO
[3625]	992	DO jj=1,ny_global
	993	DO ji=1,nx_global
	994	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
	995	ENDDO
	996	ENDDO
[3193]	997	ENDIF
[2874]	998
	999	#endif
	1000
[3193]	1001	SELECT CASE ( cd_type )
	1002	CASE ( 'T' )
	1003	grid_loc=field_loc_center
	1004	CASE ( 'F' )
	1005	grid_loc=field_loc_NEcorner
	1006	END SELECT
[2874]	1007
[3193]	1008	SELECT CASE ( NINT(psgn) )
	1009	CASE ( -1 )
	1010	field_type=field_type_vector
	1011	CASE ( 1 )
	1012	field_type=field_type_scalar
	1013	END SELECT
[2874]	1014
	1015	#if defined key_nemocice_decomp
[3193]	1016	! Ensure CICE halos are up to date
	1017	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1018	#else
[3193]	1019	! D. Scatter pcg to CICE blocks (pc) + update halos
	1020	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
[2874]	1021	#endif
	1022
	1023	END SUBROUTINE nemo2cice
	1024
	1025	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
	1026	!!---------------------------------------------------------------------
	1027	!! * ROUTINE cice2nemo *
	1028	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
	1029	#if defined key_nemocice_decomp
	1030	!!
	1031	!! NEMO and CICE PE sub domains are identical, hence
	1032	!! there is no need to gather or scatter data from
	1033	!! one PE configuration to another.
	1034	#else
	1035	!! Automatically deal with scatter/gather between
	1036	!! different processors and blocks
	1037	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
	1038	!! B. Map pcg into NEMO global array (png)
	1039	!! C. Scatter png into NEMO field (pn) for each processor
	1040	!! D. Ensure all haloes are filled in pn
	1041	#endif
	1042	!!---------------------------------------------------------------------
	1043
[3193]	1044	CHARACTER(len=1), INTENT( in ) :: &
	1045	cd_type ! nature of pn grid-point
	1046	! ! = T or F gridpoints
	1047	REAL(wp), INTENT( in ) :: &
	1048	psgn ! control of the sign change
	1049	! ! =-1 , the sign is modified following the type of b.c. used
	1050	! ! = 1 , no sign change
	1051	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	1052
	1053	#if defined key_nemocice_decomp
[3193]	1054	INTEGER (int_kind) :: &
	1055	field_type, & ! id for type of field (scalar, vector, angle)
	1056	grid_loc ! id for location on horizontal grid
	1057	! (center, NEcorner, Nface, Eface)
[2874]	1058	#else
[3193]	1059	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	1060	#endif
	1061
[3193]	1062	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
[2874]	1063
[3193]	1064	INTEGER :: ji, jj, jn ! dummy loop indices
[2874]	1065
	1066
	1067	#if defined key_nemocice_decomp
	1068
[3193]	1069	SELECT CASE ( cd_type )
	1070	CASE ( 'T' )
	1071	grid_loc=field_loc_center
	1072	CASE ( 'F' )
	1073	grid_loc=field_loc_NEcorner
	1074	END SELECT
[2874]	1075
[3193]	1076	SELECT CASE ( NINT(psgn) )
	1077	CASE ( -1 )
	1078	field_type=field_type_vector
	1079	CASE ( 1 )
	1080	field_type=field_type_scalar
	1081	END SELECT
[2874]	1082
[3193]	1083	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1084
	1085
[3193]	1086	pn(:,:)=0.0
	1087	DO jj=1,jpjm1
	1088	DO ji=1,jpim1
[3625]	1089	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
[3193]	1090	ENDDO
	1091	ENDDO
[2874]	1092
	1093	#else
	1094
[3193]	1095	! A. Gather CICE blocks (pc) into global array (pcg)
[2874]	1096
[3193]	1097	CALL gather_global(pcg, pc, 0, distrb_info)
[2874]	1098
	1099	! B. Map pcg into NEMO global array (png)
	1100
	1101	! Need to make sure this is robust to changes in NEMO halo rows....
	1102	! (may be OK but not spent much time thinking about it)
[3625]	1103	! Note that non-existent pcg elements may be used below, but
	1104	! the lbclnk call on pn will replace these with sensible values
[2874]	1105
[3193]	1106	IF (nproc==0) THEN
	1107	png(:,:,:)=0.0
	1108	DO jn=1,jpnij
[3625]	1109	DO jj=nldjt(jn),nlejt(jn)
	1110	DO ji=nldit(jn),nleit(jn)
	1111	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
[3193]	1112	ENDDO
	1113	ENDDO
	1114	ENDDO
	1115	ENDIF
[2874]	1116
[3193]	1117	! C. Scatter png into NEMO field (pn) for each processor
[2874]	1118
[3193]	1119	IF ( jpnij > 1) THEN
	1120	CALL mppsync
	1121	CALL mppscatter (png,0,pn)
	1122	CALL mppsync
	1123	ELSE
	1124	pn(:,:)=png(:,:,1)
	1125	ENDIF
[2874]	1126
	1127	#endif
	1128
[3193]	1129	! D. Ensure all haloes are filled in pn
[2874]	1130
[3193]	1131	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	1132
	1133	END SUBROUTINE cice2nemo
	1134
	1135	#else
	1136	!!----------------------------------------------------------------------
	1137	!! Default option Dummy module NO CICE sea-ice model
	1138	!!----------------------------------------------------------------------
[5215]	1139	!! $Id$
[2874]	1140	CONTAINS
	1141
[4990]	1142	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
[2874]	1143	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
	1144	END SUBROUTINE sbc_ice_cice
	1145
[4990]	1146	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
[2874]	1147	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
	1148	END SUBROUTINE cice_sbc_init
	1149
	1150	SUBROUTINE cice_sbc_final ! Dummy routine
	1151	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
	1152	END SUBROUTINE cice_sbc_final
	1153
	1154	#endif
	1155
	1156	!!======================================================================
	1157	END MODULE sbcice_cice

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

Download in other formats: