source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9631

Last change on this file since 9631 was 9631, checked in by davestorkey, 2 years ago

Bug fix in nemo2cice in sbcice_cice.F90.

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