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, 3 years ago

Working version of ice thickness assimilation updates

File size: 48.8 KB
Line 
1MODULE 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$
112CONTAINS
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_sfac*sprecip(:,:),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(:,:)=rau0*rcp*fse3t_m(:,:)*(sstfrz(:,:)-sst_m(:,:))/(2.0*dt) 
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
663taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1**2. + ztmp2**2.) 
664
665! Freshwater fluxes
666
667      IF (ksbc == jp_flx) THEN
668! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
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)
929901      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 )
933902      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$
1214CONTAINS
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   !!======================================================================
1231END MODULE sbcice_cice
Note: See TracBrowser for help on using the repository browser.