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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 10181

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

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

Download in other formats: