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

Line
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
16	USE domvvl
17	USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater
18	USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0
19	USE in_out_manager ! I/O manager
20	USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit
21	USE lib_mpp ! distributed memory computing library
22	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
23	USE wrk_nemo ! work arrays
24	USE timing ! Timing
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
40	# if defined key_cice4
41	USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
42	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
43	strocnxT,strocnyT, &
44	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, &
45	fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, &
46	flatn_f,fsurfn_f,fcondtopn_f, &
47	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
48	swvdr,swvdf,swidr,swidf,Tf
49	USE ice_therm_vertical, only: calc_Tsfc
50	#else
51	USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,&
52	vsnon,vice,vicen,nt_Tsfc
53	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
54	strocnxT,strocnyT, &
55	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai, &
56	fresh_ai,fhocn_ai,fswthru_ai,frzmlt, &
57	flatn_f,fsurfn_f,fcondtopn_f, &
58	daice_da,fresh_da,fsalt_da, &
59	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
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
69	#endif
70	USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
71	USE ice_atmo, only: calc_strair
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
85	INTEGER :: ji_off
86	INTEGER :: jj_off
87
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
109	!! $Id$
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
121	SUBROUTINE sbc_ice_cice( kt, ksbc )
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	!!
128	!! ** Method : - Get any extra forcing fields for CICE
129	!! - Prepare forcing fields
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:
136	!! utau, vtau, qns , qsr, emp , sfx
137	!!---------------------------------------------------------------------
138	INTEGER, INTENT(in) :: kt ! ocean time step
139	INTEGER, INTENT(in) :: ksbc ! surface forcing type
140	!!----------------------------------------------------------------------
141	!
142	IF( nn_timing == 1 ) CALL timing_start('sbc_ice_cice')
143	!
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
149	IF ( ksbc == jp_flx ) THEN
150	CALL cice_sbc_force(kt)
151	ELSE IF ( ksbc == jp_purecpl ) THEN
152	CALL sbc_cpl_ice_flx( 1.0-fr_i )
153	ENDIF
154
155	CALL cice_sbc_in ( kt, ksbc )
156	CALL CICE_Run
157	CALL cice_sbc_out ( kt, ksbc )
158
159	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(kt+1)
160
161	ENDIF ! End sea-ice time step only
162	!
163	IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_cice')
164
165	END SUBROUTINE sbc_ice_cice
166
167	SUBROUTINE cice_sbc_init (ksbc)
168	!!---------------------------------------------------------------------
169	!! * ROUTINE cice_sbc_init *
170	!! ** Purpose: Initialise ice related fields for NEMO and coupling
171	!!
172	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
173	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
174	REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d
175	INTEGER :: ji, jj, jl, jk ! dummy loop indices
176	!!---------------------------------------------------------------------
177
178	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_init')
179	!
180	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
181	!
182	IF(lwp) WRITE(numout,*)'cice_sbc_init'
183
184	ji_off = INT ( (jpiglo - nx_global) / 2 )
185	jj_off = INT ( (jpjglo - ny_global) / 2 )
186
187	! Initialize CICE
188	CALL CICE_Initialize
189
190	! Do some CICE consistency checks
191	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
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
195	ELSEIF (ksbc == jp_core) THEN
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
200
201
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' )
204	IF( sbc_ice_cice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
205
206	! Ensure that no temperature points are below freezing if not a NEMO restart
207	IF( .NOT. ln_rstart ) THEN
208
209	CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d )
210	DO jk=1,jpk
211	CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), fsdept_n(:,:,jk) )
212	ENDDO
213	tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d )
214	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
215	CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d )
216
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
223
224	ENDIF
225
226	! calculate surface freezing temperature and send to CICE
227	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1))
228	CALL nemo2cice(sstfrz,Tf, 'T', 1. )
229
230	CALL cice2nemo(aice,fr_i, 'T', 1. )
231	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
232	DO jl=1,ncat
233	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
234	ENDDO
235	ENDIF
236
237	! T point to U point
238	! T point to V point
239	fr_iu(:,:)=0.0
240	fr_iv(:,:)=0.0
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
247
248	CALL lbc_lnk ( fr_iu , 'U', 1. )
249	CALL lbc_lnk ( fr_iv , 'V', 1. )
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
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
301	ENDIF
302
303	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
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
310	!
311	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_init')
312	!
313	END SUBROUTINE cice_sbc_init
314
315
316	SUBROUTINE cice_sbc_in (kt, ksbc)
317	!!---------------------------------------------------------------------
318	!! * ROUTINE cice_sbc_in *
319	!! ** Purpose: Set coupling fields and pass to CICE
320	!!---------------------------------------------------------------------
321	INTEGER, INTENT(in ) :: kt ! ocean time step
322	INTEGER, INTENT(in ) :: ksbc ! surface forcing type
323
324	INTEGER :: ji, jj, jl ! dummy loop indices
325	REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice
326	REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
327	REAL(wp) :: zintb, zintn ! dummy argument
328	!!---------------------------------------------------------------------
329
330	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_in')
331	!
332	CALL wrk_alloc( jpi,jpj, ztmp, zpice )
333	CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
334
335	IF( kt == nit000 ) THEN
336	IF(lwp) WRITE(numout,*)'cice_sbc_in'
337	ENDIF
338
339	ztmp(:,:)=0.0
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
346	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
347
348	ztmpn(:,:,:)=0.0
349
350	! x comp of wind stress (CI_1)
351	! U point to F point
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. )
359
360	! y comp of wind stress (CI_2)
361	! V point to F point
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. )
369
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
382	! Surface downward latent heat flux (CI_5)
383	IF (ksbc == jp_flx) THEN
384	DO jl=1,ncat
385	ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
386	ENDDO
387	ELSE IF (ksbc == jp_purecpl) THEN
388	DO jl=1,ncat
389	ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl)
390	ENDDO
391	ELSE
392	!In coupled mode - qla_ice calculated in sbc_cpl for each category
393	ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat)
394	ENDIF
395
396	DO jl=1,ncat
397	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
398
399	! GBM conductive flux through ice (CI_6)
400	! Convert to GBM
401	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
402	ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
403	ELSE
404	ztmp(:,:) = botmelt(:,:,jl)
405	ENDIF
406	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
407
408	! GBM surface heat flux (CI_7)
409	! Convert to GBM
410	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
411	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)
412	ELSE
413	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
414	ENDIF
415	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
416	ENDDO
417
418	ELSE IF (ksbc == jp_core) THEN
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)
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)
435	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
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)
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)
448	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
449	CALL nemo2cice(ztmp,swvdr,'T', 1. )
450	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
451	CALL nemo2cice(ztmp,swvdf,'T', 1. )
452	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
453	CALL nemo2cice(ztmp,swidr,'T', 1. )
454	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
455	CALL nemo2cice(ztmp,swidf,'T', 1. )
456
457	ENDIF
458
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
465	! Snowfall
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
468	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
469	CALL nemo2cice(ztmp,fsnow,'T', 1. )
470
471	! Rainfall
472	IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
473	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
474	CALL nemo2cice(ztmp,frain,'T', 1. )
475
476	! Recalculate freezing temperature and send to CICE
477	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1))
478	CALL nemo2cice(sstfrz,Tf,'T', 1. )
479
480	! Freezing/melting potential
481	! Calculated over NEMO leapfrog timestep (hence 2*dt)
482	nfrzmlt(:,:)=rau0rcpfse3t_m(:,:)(sstfrz(:,:)-sst_m(:,:))/(2.0dt)
483	CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. )
484
485	! SST and SSS
486
487	CALL nemo2cice(sst_m,sst,'T', 1. )
488	CALL nemo2cice(sss_m,sss,'T', 1. )
489
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
497	! x comp and y comp of surface ocean current
498	! U point to F point
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. )
505
506	! V point to F point
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. )
513
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
531	! x comp and y comp of sea surface slope (on F points)
532	! T point to F point
533	DO jj=1,jpjm1
534	DO ji=1,jpim1
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) ) &
537	* fmask(ji,jj,1)
538	ENDDO
539	ENDDO
540	CALL nemo2cice(ztmp,ss_tltx,'F', -1. )
541
542	! T point to F point
543	DO jj=1,jpjm1
544	DO ji=1,jpim1
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) ) &
547	* fmask(ji,jj,1)
548	ENDDO
549	ENDDO
550	CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
551
552	CALL wrk_dealloc( jpi,jpj, ztmp, zpice )
553	CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
554	!
555	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_in')
556	!
557	END SUBROUTINE cice_sbc_in
558
559
560	SUBROUTINE cice_sbc_out (kt,ksbc)
561	!!---------------------------------------------------------------------
562	!! * ROUTINE cice_sbc_out *
563	!! ** Purpose: Get fields from CICE and set surface fields for NEMO
564	!!---------------------------------------------------------------------
565	INTEGER, INTENT( in ) :: kt ! ocean time step
566	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
567
568	INTEGER :: ji, jj, jl ! dummy loop indices
569	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
570	!!---------------------------------------------------------------------
571
572	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_out')
573	!
574	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
575
576	IF( kt == nit000 ) THEN
577	IF(lwp) WRITE(numout,*)'cice_sbc_out'
578	ENDIF
579
580	! x comp of ocean-ice stress
581	CALL cice2nemo(strocnx,ztmp1,'F', -1. )
582	ss_iou(:,:)=0.0
583	! F point to U point
584	DO jj=2,jpjm1
585	DO ji=2,jpim1
586	ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
587	ENDDO
588	ENDDO
589	CALL lbc_lnk( ss_iou , 'U', -1. )
590
591	! y comp of ocean-ice stress
592	CALL cice2nemo(strocny,ztmp1,'F', -1. )
593	ss_iov(:,:)=0.0
594	! F point to V point
595
596	DO jj=1,jpjm1
597	DO ji=2,jpim1
598	ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
599	ENDDO
600	ENDDO
601	CALL lbc_lnk( ss_iov , 'V', -1. )
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]
606	! strocnx and strocny already weighted by ice fraction in CICE so not done here
607
608	utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
609	vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
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.)
619
620	! Freshwater fluxes
621
622	IF (ksbc == jp_flx) THEN
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)
627	emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
628	ELSE IF (ksbc == jp_core) THEN
629	emp(:,:) = (1.0-fr_i(:,:))*emp(:,:)
630	ELSE IF (ksbc == jp_purecpl) THEN
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
633	emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
634	ENDIF
635
636	#if defined key_cice4
637	CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
638	CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
639	#else
640	CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
641	CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
642	#endif
643
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(:,:)
652	fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
653
654	CALL lbc_lnk( emp , 'T', 1. )
655	CALL lbc_lnk( sfx , 'T', 1. )
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
661	IF (ksbc == jp_core) THEN
662	qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
663	qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
664	ENDIF
665	! Take into account snow melting except for fully coupled when already in qns_tot
666	IF (ksbc == jp_purecpl) THEN
667	qsr(:,:)= qsr_tot(:,:)
668	qns(:,:)= qns_tot(:,:)
669	ELSE
670	qns(:,:)= qns(:,:)-sprecip(:,:)Lfresh(1.0-fr_i(:,:))
671	ENDIF
672
673	! Now add in ice / snow related terms
674	! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
675	#if defined key_cice4
676	CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
677	#else
678	CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
679	#endif
680	qsr(:,:)=qsr(:,:)+ztmp1(:,:)
681	CALL lbc_lnk( qsr , 'T', 1. )
682
683	DO jj=1,jpj
684	DO ji=1,jpi
685	nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
686	ENDDO
687	ENDDO
688
689	#if defined key_cice4
690	CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
691	#else
692	CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
693	#endif
694	qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
695
696	CALL lbc_lnk( qns , 'T', 1. )
697
698	! Prepare for the following CICE time-step
699
700	CALL cice2nemo(aice,fr_i,'T', 1. )
701	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
702	DO jl=1,ncat
703	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
704	ENDDO
705	ENDIF
706
707	! T point to U point
708	! T point to V point
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
715
716	CALL lbc_lnk ( fr_iu , 'U', 1. )
717	CALL lbc_lnk ( fr_iv , 'V', 1. )
718
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
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
734	! Release work space
735
736	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
737	!
738	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_out')
739	!
740	END SUBROUTINE cice_sbc_out
741
742
743	SUBROUTINE cice_sbc_hadgam( kt )
744	!!---------------------------------------------------------------------
745	!! * ROUTINE cice_sbc_hadgam *
746	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
747	!!
748	!!
749	INTEGER, INTENT( in ) :: kt ! ocean time step
750	!!---------------------------------------------------------------------
751
752	INTEGER :: jl ! dummy loop index
753	INTEGER :: ierror
754
755	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_hadgam')
756	!
757	IF( kt == nit000 ) THEN
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' )
760	ENDIF
761
762	! ! =========================== !
763	! ! Prepare Coupling fields !
764	! ! =========================== !
765
766	! x and y comp of ice velocity
767
768	CALL cice2nemo(uvel,u_ice,'F', -1. )
769	CALL cice2nemo(vvel,v_ice,'F', -1. )
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
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. )
778	ENDDO
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
807	!
808	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_hadgam')
809	!
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
821	CALL CICE_Finalize
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
860	INTEGER :: ios
861	!!---------------------------------------------------------------------
862
863	! ! ====================== !
864	IF( kt == nit000 ) THEN ! First call kt=nit000 !
865	! ! ====================== !
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
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 )
890
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 )
894	IF(lwm) WRITE ( numond, namsbc_cice )
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)
925	! May be better to do this conversion somewhere else
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
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
974	#if !defined key_nemocice_decomp
975	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
976	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
977	#endif
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
982	! (center, NEcorner, Nface, Eface)
983
984	INTEGER :: ji, jj, jn ! dummy loop indices
985
986	! A. Ensure all haloes are filled in NEMO field (pn)
987
988	CALL lbc_lnk( pn , cd_type, psgn )
989
990	#if defined key_nemocice_decomp
991
992	! Copy local domain data from NEMO to CICE field
993	pc(:,:,1)=0.0
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)
997	ENDDO
998	ENDDO
999
1000	#else
1001
1002	! B. Gather pn into global array (png)
1003
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
1011
1012	! C. Map png into CICE global array (pcg)
1013
1014	! Need to make sure this is robust to changes in NEMO halo rows....
1015	! (may be OK but not 100% sure)
1016
1017	IF (nproc==0) THEN
1018	! pcg(:,:)=0.0
1019	DO jn=1,jpnij
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)
1023	ENDDO
1024	ENDDO
1025	ENDDO
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
1031	ENDIF
1032
1033	#endif
1034
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
1041
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
1048
1049	#if defined key_nemocice_decomp
1050	! Ensure CICE halos are up to date
1051	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1052	#else
1053	! D. Scatter pcg to CICE blocks (pc) + update halos
1054	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
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
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
1086
1087	#if defined key_nemocice_decomp
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)
1092	#else
1093	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
1094	#endif
1095
1096	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
1097
1098	INTEGER :: ji, jj, jn ! dummy loop indices
1099
1100
1101	#if defined key_nemocice_decomp
1102
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
1109
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
1116
1117	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1118
1119
1120	pn(:,:)=0.0
1121	DO jj=1,jpjm1
1122	DO ji=1,jpim1
1123	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
1124	ENDDO
1125	ENDDO
1126
1127	#else
1128
1129	! A. Gather CICE blocks (pc) into global array (pcg)
1130
1131	CALL gather_global(pcg, pc, 0, distrb_info)
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)
1137	! Note that non-existent pcg elements may be used below, but
1138	! the lbclnk call on pn will replace these with sensible values
1139
1140	IF (nproc==0) THEN
1141	png(:,:,:)=0.0
1142	DO jn=1,jpnij
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)
1146	ENDDO
1147	ENDDO
1148	ENDDO
1149	ENDIF
1150
1151	! C. Scatter png into NEMO field (pn) for each processor
1152
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
1160
1161	#endif
1162
1163	! D. Ensure all haloes are filled in pn
1164
1165	CALL lbc_lnk( pn , cd_type, psgn )
1166
1167	END SUBROUTINE cice2nemo
1168
1169	#else
1170	!!----------------------------------------------------------------------
1171	!! Default option Dummy module NO CICE sea-ice model
1172	!!----------------------------------------------------------------------
1173	!! $Id$
1174	CONTAINS
1175
1176	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
1177	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
1178	END SUBROUTINE sbc_ice_cice
1179
1180	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
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: