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

Last change on this file since 9819 was 9819, checked in by dancopsey, 2 years ago

Added a bunch of fixes and a nemo2cice print statement.

File size: 46.9 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.,'sst1') 
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.,'tf1' )
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      WRITE(numout,*)'zevap_ice 7: min, max = ', MINVAL(zevap_ice(:,:,:)), MAXVAL(zevap_ice(:,:,:))
323
324      ztmp(:,:)=0.0
325
326! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
327! the first time-step)
328
329! forced and coupled case
330
331      IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
332
333         ztmpn(:,:,:)=0.0
334
335! x comp of wind stress (CI_1)
336! U point to F point
337         DO jj=1,jpjm1
338            DO ji=1,jpi
339               ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
340                                    + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
341            ENDDO
342         ENDDO
343         
344         CALL nemo2cice(ztmp,strax,'F', -1.,'strax' )         
345
346! y comp of wind stress (CI_2)
347! V point to F point
348         DO jj=1,jpj
349            DO ji=1,jpim1
350               ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
351                                    + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
352            ENDDO
353         ENDDO
354         CALL nemo2cice(ztmp,stray,'F', -1.,'stray' )
355
356
357! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in
358! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby
359! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing
360! gridbox mean fluxes in the UM by future ice concentration obtained through 
361! OASIS.  This allows for a much more realistic apportionment of energy through
362! the ice - and conserves energy.
363! Therefore the fluxes are now divided by ice concentration in the coupled
364! formulation (jp_purecpl) as well as for jp_flx.  This NEMO branch should only
365! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at
366! which point the GSI8 UM changes were committed.
367
368! Surface downward latent heat flux (CI_5)
369         IF (ksbc == jp_flx) THEN
370            DO jl=1,ncat
371               ztmpn(:,:,jl)=zevap_ice(:,:,1)*a_i(:,:,jl) * (-1.0*Lsub)
372            ENDDO
373         ELSE IF (ksbc == jp_purecpl) THEN
374            DO jl=1,ncat
375               ztmpn(:,:,jl)=zevap_ice(:,:,jl)*a_i(:,:,jl) * (-1.0*Lsub)
376            ENDDO
377    ELSE
378           !In coupled mode - zevap_ice calculated in sbc_cpl for each category
379           ztmpn(:,:,1:ncat)=zevap_ice(:,:,1:ncat) * (-1.0*Lsub)
380         ENDIF
381
382         DO jl=1,ncat
383            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1.,'flatn' )
384
385! GBM conductive flux through ice (CI_6)
386!  Convert to GBM
387            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
388               ztmp(:,:) = qcn_ice(:,:,jl)*a_i(:,:,jl)
389            ELSE
390               ztmp(:,:) = qcn_ice(:,:,jl)
391            ENDIF
392            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1.,'fcondtopn' )
393
394! GBM surface heat flux (CI_7)
395!  Convert to GBM
396            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
397               ztmp(:,:) = (qml_ice(:,:,jl)+qcn_ice(:,:,jl))*a_i(:,:,jl) 
398            ELSE
399               ztmp(:,:) = (qml_ice(:,:,jl)+qcn_ice(:,:,jl))
400            ENDIF
401            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. ,'fsurfn_f')
402         ENDDO
403
404      ELSE IF (ksbc == jp_blk) THEN
405
406! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
407! x comp and y comp of atmosphere surface wind (CICE expects on T points)
408         ztmp(:,:) = wndi_ice(:,:)
409         CALL nemo2cice(ztmp,uatm,'T', -1. ,'uatm')
410         ztmp(:,:) = wndj_ice(:,:)
411         CALL nemo2cice(ztmp,vatm,'T', -1.,'vatm' )
412         ztmp(:,:) = SQRT ( wndi_ice(:,:)**2 + wndj_ice(:,:)**2 )
413         CALL nemo2cice(ztmp,wind,'T', 1.,'wind' )    ! Wind speed (m/s)
414         ztmp(:,:) = qsr_ice(:,:,1)
415         CALL nemo2cice(ztmp,fsw,'T', 1.,'fsw' )     ! Incoming short-wave (W/m^2)
416         ztmp(:,:) = qlw_ice(:,:,1)
417         CALL nemo2cice(ztmp,flw,'T', 1.,'flw' )     ! Incoming long-wave (W/m^2)
418         ztmp(:,:) = tatm_ice(:,:)
419         CALL nemo2cice(ztmp,Tair,'T', 1.,'tair' )    ! Air temperature (K)
420         CALL nemo2cice(ztmp,potT,'T', 1.,'pott' )    ! Potential temp (K)
421! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows 
422         ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )   
423                                               ! Constant (101000.) atm pressure assumed
424         CALL nemo2cice(ztmp,rhoa,'T', 1.,'rhoa' )    ! Air density (kg/m^3)
425         ztmp(:,:) = qatm_ice(:,:)
426         CALL nemo2cice(ztmp,Qa,'T', 1.,'qa' )      ! Specific humidity (kg/kg)
427         ztmp(:,:)=10.0
428         CALL nemo2cice(ztmp,zlvl,'T', 1.,'zlvl' )    ! Atmos level height (m)
429
430! May want to check all values are physically realistic (as in CICE routine
431! prepare_forcing)?
432
433! Divide shortwave into spectral bands (as in prepare_forcing)
434         ztmp(:,:)=qsr_ice(:,:,1)*frcvdr       ! visible direct
435         CALL nemo2cice(ztmp,swvdr,'T', 1.,'swvdr' )             
436         ztmp(:,:)=qsr_ice(:,:,1)*frcvdf       ! visible diffuse
437         CALL nemo2cice(ztmp,swvdf,'T', 1.,'swvdf' )             
438         ztmp(:,:)=qsr_ice(:,:,1)*frcidr       ! near IR direct
439         CALL nemo2cice(ztmp,swidr,'T', 1.,'swidr' )
440         ztmp(:,:)=qsr_ice(:,:,1)*frcidf       ! near IR diffuse
441         CALL nemo2cice(ztmp,swidf,'T', 1. ,'swidf')
442
443      ENDIF
444
445#if defined key_asminc
446!Ice concentration change (from assimilation)
447      ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
448      Call nemo2cice(ztmp,daice_da,'T', 1. ,'daice_da')
449#endif 
450
451! Snowfall
452! Ensure fsnow is positive (as in CICE routine prepare_forcing)
453      IF( iom_use('snowpre') )   CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit 
454      ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0) 
455      CALL nemo2cice(ztmp,fsnow,'T', 1.,'fsnow' ) 
456
457! Rainfall
458      IF( iom_use('precip') )   CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
459      ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
460      CALL nemo2cice(ztmp,frain,'T', 1.,'frain' ) 
461
462! Recalculate freezing temperature and send to CICE
463      CALL eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1)) 
464      CALL nemo2cice(sstfrz,Tf,'T', 1. ,'Tf')
465
466! Freezing/melting potential
467! Calculated over NEMO leapfrog timestep (hence 2*dt)
468      nfrzmlt(:,:)=rau0*rcp*e3t_m(:,:)*(sstfrz(:,:)-sst_m(:,:))/(2.0*dt) 
469      CALL nemo2cice(nfrzmlt,frzmlt,'T', 1.,'frzmlt' )
470
471! SST  and SSS
472
473      CALL nemo2cice(sst_m,sst,'T', 1.,'sst' )
474      CALL nemo2cice(sss_m,sss,'T', 1.,'sss' )
475
476      IF( ksbc == jp_purecpl ) THEN
477! Sea ice surface skin temperature
478         DO jl=1,ncat
479           CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1., 'trcrn')
480         ENDDO 
481      ENDIF
482
483! x comp and y comp of surface ocean current
484! U point to F point
485      DO jj=1,jpjm1
486         DO ji=1,jpi
487            ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
488         ENDDO
489      ENDDO
490      CALL nemo2cice(ztmp,uocn,'F', -1., 'uocn' )
491
492! V point to F point
493      DO jj=1,jpj
494         DO ji=1,jpim1
495            ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
496         ENDDO
497      ENDDO
498      CALL nemo2cice(ztmp,vocn,'F', -1. ,'vocn')
499
500      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
501          !
502          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
503          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
504         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
505          !
506          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
507          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
508         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
509          !
510         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
511          !
512         !
513      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
514         zpice(:,:) = ssh_m(:,:)
515      ENDIF
516
517! x comp and y comp of sea surface slope (on F points)
518! T point to F point
519      DO jj = 1, jpjm1
520         DO ji = 1, jpim1
521            ztmp(ji,jj)=0.5 * (  (zpice(ji+1,jj  )-zpice(ji,jj  )) * r1_e1u(ji,jj  )    &
522               &               + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1)  ) * fmask(ji,jj,1)
523         END DO
524      END DO
525      CALL nemo2cice( ztmp,ss_tltx,'F', -1. , 'ss_tltx')
526
527! T point to F point
528      DO jj = 1, jpjm1
529         DO ji = 1, jpim1
530            ztmp(ji,jj)=0.5 * (  (zpice(ji  ,jj+1)-zpice(ji  ,jj)) * r1_e2v(ji  ,jj)    &
531               &               + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj)  ) *  fmask(ji,jj,1)
532         END DO
533      END DO
534      CALL nemo2cice(ztmp,ss_tlty,'F', -1. ,'ss_tlty')
535      !
536   END SUBROUTINE cice_sbc_in
537
538
539   SUBROUTINE cice_sbc_out( kt, ksbc )
540      !!---------------------------------------------------------------------
541      !!                    ***  ROUTINE cice_sbc_out  ***
542      !! ** Purpose: Get fields from CICE and set surface fields for NEMO
543      !!---------------------------------------------------------------------
544      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
545      INTEGER, INTENT( in  ) ::   ksbc ! surface forcing type
546     
547      INTEGER  ::   ji, jj, jl                 ! dummy loop indices
548      REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
549      !!---------------------------------------------------------------------
550      !
551      IF( kt == nit000 )  THEN
552         IF(lwp) WRITE(numout,*)'cice_sbc_out'
553      ENDIF
554     
555! x comp of ocean-ice stress
556      CALL cice2nemo(strocnx,ztmp1,'F', -1. )
557      ss_iou(:,:)=0.0
558! F point to U point
559      DO jj=2,jpjm1
560         DO ji=2,jpim1
561            ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
562         ENDDO
563      ENDDO
564      CALL lbc_lnk( ss_iou , 'U', -1. )
565
566! y comp of ocean-ice stress
567      CALL cice2nemo(strocny,ztmp1,'F', -1. )
568      ss_iov(:,:)=0.0
569! F point to V point
570
571      DO jj=1,jpjm1
572         DO ji=2,jpim1
573            ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
574         ENDDO
575      ENDDO
576      CALL lbc_lnk( ss_iov , 'V', -1. )
577
578! x and y comps of surface stress
579! Combine wind stress and ocean-ice stress
580! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
581! strocnx and strocny already weighted by ice fraction in CICE so not done here
582
583      utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
584      vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)     
585 
586! Also need ice/ocean stress on T points so that taum can be updated
587! This interpolation is already done in CICE so best to use those values
588      CALL cice2nemo(strocnxT,ztmp1,'T',-1.) 
589      CALL cice2nemo(strocnyT,ztmp2,'T',-1.) 
590 
591! Update taum with modulus of ice-ocean stress
592! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
593taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2) 
594
595! Freshwater fluxes
596
597      IF (ksbc == jp_flx) THEN
598! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
599! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below
600! Not ideal since aice won't be the same as in the atmosphere. 
601! Better to use evap and tprecip? (but for now don't read in evap in this case)
602         emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
603      ELSE IF (ksbc == jp_blk) THEN
604         emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)       
605      ELSE IF (ksbc == jp_purecpl) THEN
606! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
607! This is currently as required with the coupling fields from the UM atmosphere
608         emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:) 
609      ENDIF
610
611#if defined key_cice4
612      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
613      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
614#else
615      CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
616      CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
617#endif
618
619! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
620! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
621! which has to be compensated for by the ocean salinity potentially going negative
622! This check breaks conservation but seems reasonable until we have prognostic ice salinity
623! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
624      WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
625      sfx(:,:)=ztmp2(:,:)*1000.0
626      emp(:,:)=emp(:,:)-ztmp1(:,:)
627      fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
628     
629      CALL lbc_lnk_multi( emp , 'T', 1., sfx , 'T', 1. )
630
631! Solar penetrative radiation and non solar surface heat flux
632
633! Scale qsr and qns according to ice fraction (bulk formulae only)
634
635      IF (ksbc == jp_blk) THEN
636         qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
637         qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
638      ENDIF
639! Take into account snow melting except for fully coupled when already in qns_tot
640      IF (ksbc == jp_purecpl) THEN
641         qsr(:,:)= qsr_tot(:,:)
642         qns(:,:)= qns_tot(:,:)
643      ELSE
644         qns(:,:)= qns(:,:)-sprecip(:,:)*Lfresh*(1.0-fr_i(:,:))
645      ENDIF
646
647! Now add in ice / snow related terms
648! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
649#if defined key_cice4
650      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
651#else
652      CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
653#endif
654      qsr(:,:)=qsr(:,:)+ztmp1(:,:)
655      CALL lbc_lnk( qsr , 'T', 1. )
656
657      DO jj=1,jpj
658         DO ji=1,jpi
659            nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
660         ENDDO
661      ENDDO
662
663#if defined key_cice4
664      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
665#else
666      CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
667#endif
668      qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
669
670      CALL lbc_lnk( qns , 'T', 1. )
671
672! Prepare for the following CICE time-step
673
674      CALL cice2nemo(aice,fr_i,'T', 1. )
675      IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
676         DO jl=1,ncat
677            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
678         ENDDO
679      ENDIF
680
681! T point to U point
682! T point to V point
683      DO jj=1,jpjm1
684         DO ji=1,jpim1
685            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
686            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
687         ENDDO
688      ENDDO
689
690      CALL lbc_lnk_multi( fr_iu , 'U', 1., fr_iv , 'V', 1. )
691
692      ! set the snow+ice mass
693      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
694      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
695      snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
696      snwice_mass_b(:,:) = snwice_mass(:,:)
697      snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
698      !
699   END SUBROUTINE cice_sbc_out
700
701
702   SUBROUTINE cice_sbc_hadgam( kt )
703      !!---------------------------------------------------------------------
704      !!                    ***  ROUTINE cice_sbc_hadgam  ***
705      !! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
706      !!
707      !!
708      !!---------------------------------------------------------------------
709      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
710      !!
711      INTEGER  ::   jl                        ! dummy loop index
712      INTEGER  ::   ierror
713      !!---------------------------------------------------------------------
714#if defined key_asminc
715! Import fresh water and salt flux due to seaice da
716      CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
717      CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
718#endif
719
720      !
721      !                                         ! =========================== !
722      !                                         !   Prepare Coupling fields   !
723      !                                         ! =========================== !
724      !
725      ! x and y comp of ice velocity
726      !
727      CALL cice2nemo(uvel,u_ice,'F', -1. )
728      CALL cice2nemo(vvel,v_ice,'F', -1. )
729      !
730      ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out 
731      !
732      ! Snow and ice thicknesses (CO_2 and CO_3)
733      !
734      DO jl = 1, ncat
735         CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
736         CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
737      END DO
738
739#if ! defined key_cice4
740! Meltpond fraction and depth
741      DO jl = 1,ncat
742         CALL cice2nemo(apeffn(:,:,jl,:),a_ip(:,:,jl),'T', 1. )
743         CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),v_ip(:,:,jl),'T', 1. )
744      ENDDO
745     
746      v_ip(:,:,:) = v_ip(:,:,:) * a_i(:,:,:)
747#endif
748
749
750! If using multilayers thermodynamics in CICE then get top layer temperature
751! and effective conductivity       
752!! When using NEMO with CICE, this change requires use of
753!! one of the following two CICE branches:
754!! - at CICE5.0,   hadax/r1015_GSI8_with_GSI7
755!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
756      IF (heat_capacity) THEN
757         DO jl = 1,ncat
758            CALL cice2nemo(Tn_top(:,:,jl,:),t1_ice(:,:,jl),'T', 1. )
759            CALL cice2nemo(keffn_top(:,:,jl,:),cnd_ice(:,:,jl),'T', 1. )
760         ENDDO
761! Convert surface temperature to Kelvin
762         t1_ice(:,:,:)=t1_ice(:,:,:)+rt0
763      ELSE
764         t1_ice(:,:,:) = 0.0
765         cnd_ice(:,:,:) = 0.0
766      ENDIF       
767
768      !
769   END SUBROUTINE cice_sbc_hadgam
770
771
772   SUBROUTINE cice_sbc_final
773      !!---------------------------------------------------------------------
774      !!                    ***  ROUTINE cice_sbc_final  ***
775      !! ** Purpose: Finalize CICE
776      !!---------------------------------------------------------------------
777      !
778      IF(lwp) WRITE(numout,*)'cice_sbc_final'
779      !
780      CALL CICE_Finalize
781      !
782   END SUBROUTINE cice_sbc_final
783
784
785   SUBROUTINE cice_sbc_force (kt)
786      !!---------------------------------------------------------------------
787      !!                    ***  ROUTINE cice_sbc_force  ***
788      !! ** Purpose : Provide CICE forcing from files
789      !!
790      !!---------------------------------------------------------------------
791      !! ** Method  :   READ monthly flux file in NetCDF files
792      !!     
793      !!  snowfall   
794      !!  rainfall   
795      !!  sublimation rate   
796      !!  topmelt (category)
797      !!  botmelt (category)
798      !!
799      !! History :
800      !!----------------------------------------------------------------------
801      USE iom
802      !!
803      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
804      !!
805      INTEGER  ::   ierror             ! return error code
806      INTEGER  ::   ifpr               ! dummy loop index
807      !!
808      CHARACTER(len=100) ::  cn_dir                            !   Root directory for location of CICE forcing files
809      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
810      TYPE(FLD_N) ::   sn_snow, sn_rain, sn_sblm               ! informations about the fields to be read
811      TYPE(FLD_N) ::   sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
812      TYPE(FLD_N) ::   sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5 
813      !!
814      NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm,   &
815         &                          sn_top1, sn_top2, sn_top3, sn_top4, sn_top5,  &
816         &                          sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
817      INTEGER :: ios
818      !!---------------------------------------------------------------------
819
820      !                                         ! ====================== !
821      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
822         !                                      ! ====================== !
823         ! namsbc_cice is not yet in the reference namelist
824         ! set file information (default values)
825         cn_dir = './'       ! directory in which the model is executed
826
827         ! (NB: frequency positive => hours, negative => months)
828         !            !    file          ! frequency !  variable    ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! landmask
829         !            !    name          !  (hours)  !   name       !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! file
830         sn_snow = FLD_N( 'snowfall_1m'  ,    -1.    ,  'snowfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    ) 
831         sn_rain = FLD_N( 'rainfall_1m'  ,    -1.    ,  'rainfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    ) 
832         sn_sblm = FLD_N( 'sublim_1m'    ,    -1.    ,  'sublim'    ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
833         sn_top1 = FLD_N( 'topmeltn1_1m' ,    -1.    ,  'topmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
834         sn_top2 = FLD_N( 'topmeltn2_1m' ,    -1.    ,  'topmeltn2' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
835         sn_top3 = FLD_N( 'topmeltn3_1m' ,    -1.    ,  'topmeltn3' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
836         sn_top4 = FLD_N( 'topmeltn4_1m' ,    -1.    ,  'topmeltn4' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
837         sn_top5 = FLD_N( 'topmeltn5_1m' ,    -1.    ,  'topmeltn5' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
838         sn_bot1 = FLD_N( 'botmeltn1_1m' ,    -1.    ,  'botmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
839         sn_bot2 = FLD_N( 'botmeltn2_1m' ,    -1.    ,  'botmeltn2' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
840         sn_bot3 = FLD_N( 'botmeltn3_1m' ,    -1.    ,  'botmeltn3' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
841         sn_bot4 = FLD_N( 'botmeltn4_1m' ,    -1.    ,  'botmeltn4' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
842         sn_bot5 = FLD_N( 'botmeltn5_1m' ,    -1.    ,  'botmeltn5' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
843
844         REWIND( numnam_ref )              ! Namelist namsbc_cice in reference namelist :
845         READ  ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
846901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
847
848         REWIND( numnam_cfg )              ! Namelist namsbc_cice in configuration namelist : Parameters of the run
849         READ  ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
850902      IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
851         IF(lwm) WRITE ( numond, namsbc_cice )
852
853         ! store namelist information in an array
854         slf_i(jp_snow) = sn_snow   ;   slf_i(jp_rain) = sn_rain   ;   slf_i(jp_sblm) = sn_sblm
855         slf_i(jp_top1) = sn_top1   ;   slf_i(jp_top2) = sn_top2   ;   slf_i(jp_top3) = sn_top3
856         slf_i(jp_top4) = sn_top4   ;   slf_i(jp_top5) = sn_top5   ;   slf_i(jp_bot1) = sn_bot1
857         slf_i(jp_bot2) = sn_bot2   ;   slf_i(jp_bot3) = sn_bot3   ;   slf_i(jp_bot4) = sn_bot4
858         slf_i(jp_bot5) = sn_bot5
859         
860         ! set sf structure
861         ALLOCATE( sf(jpfld), STAT=ierror )
862         IF( ierror > 0 ) THEN
863            CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' )   ;   RETURN
864         ENDIF
865
866         DO ifpr= 1, jpfld
867            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
868            ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
869         END DO
870
871         ! fill sf with slf_i and control print
872         CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
873         !
874      ENDIF
875
876      CALL fld_read( kt, nn_fsbc, sf )           ! Read input fields and provides the
877      !                                          ! input fields at the current time-step
878
879      ! set the fluxes from read fields
880      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
881      tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
882! May be better to do this conversion somewhere else
883      zevap_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
884      qml_ice(:,:,1) = sf(jp_top1)%fnow(:,:,1)
885      qml_ice(:,:,2) = sf(jp_top2)%fnow(:,:,1)
886      qml_ice(:,:,3) = sf(jp_top3)%fnow(:,:,1)
887      qml_ice(:,:,4) = sf(jp_top4)%fnow(:,:,1)
888      qml_ice(:,:,5) = sf(jp_top5)%fnow(:,:,1)
889      qcn_ice(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
890      qcn_ice(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
891      qcn_ice(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
892      qcn_ice(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
893      qcn_ice(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
894
895      ! control print (if less than 100 time-step asked)
896      IF( nitend-nit000 <= 100 .AND. lwp ) THEN
897         WRITE(numout,*) 
898         WRITE(numout,*) '        read forcing fluxes for CICE OK'
899         CALL FLUSH(numout)
900      ENDIF
901
902   END SUBROUTINE cice_sbc_force
903
904   SUBROUTINE nemo2cice( pn, pc, cd_type, psgn, varname)
905      !!---------------------------------------------------------------------
906      !!                    ***  ROUTINE nemo2cice  ***
907      !! ** Purpose :   Transfer field in NEMO array to field in CICE array. 
908#if defined key_nemocice_decomp
909      !!             
910      !!                NEMO and CICE PE sub domains are identical, hence
911      !!                there is no need to gather or scatter data from
912      !!                one PE configuration to another.
913#else
914      !!                Automatically gather/scatter between
915      !!                different processors and blocks
916      !! ** Method :    A. Ensure all haloes are filled in NEMO field (pn)
917      !!                B. Gather pn into global array (png)
918      !!                C. Map png into CICE global array (pcg)
919      !!                D. Scatter pcg to CICE blocks (pc) + update haloes 
920#endif
921      !!---------------------------------------------------------------------
922      CHARACTER(len=1), INTENT( in ) ::   &
923          cd_type       ! nature of pn grid-point
924          !             !   = T or F gridpoints
925      REAL(wp), INTENT( in ) ::   &
926          psgn          ! control of the sign change
927          !             !   =-1 , the sign is modified following the type of b.c. used
928          !             !   = 1 , no sign change
929      REAL(wp), DIMENSION(jpi,jpj) :: pn
930#if !defined key_nemocice_decomp
931      REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
932      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
933#endif
934      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
935      CHARACTER(len=10), INTENT( in ) :: varname
936      INTEGER (int_kind) :: &
937         field_type,        &! id for type of field (scalar, vector, angle)
938         grid_loc            ! id for location on horizontal grid
939                            !  (center, NEcorner, Nface, Eface)
940
941      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
942      !!---------------------------------------------------------------------
943
944!     A. Ensure all haloes are filled in NEMO field (pn)
945
946      CALL lbc_lnk( pn , cd_type, psgn )
947
948#if defined key_nemocice_decomp
949
950      ! Copy local domain data from NEMO to CICE field
951      pc(:,:,1)=0.0
952      DO jj=2,jpj-1
953         DO ji=2,jpi-1
954            pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
955         ENDDO
956      ENDDO
957
958#else
959
960!     B. Gather pn into global array (png)
961
962      IF ( jpnij > 1) THEN
963         CALL mppsync
964         CALL mppgather (pn,0,png) 
965         CALL mppsync
966      ELSE
967         png(:,:,1)=pn(:,:)
968      ENDIF
969
970!     C. Map png into CICE global array (pcg)
971
972! Need to make sure this is robust to changes in NEMO halo rows....
973! (may be OK but not 100% sure)
974
975      IF (nproc==0) THEN     
976!        pcg(:,:)=0.0
977         DO jn=1,jpnij
978            DO jj=nldjt(jn),nlejt(jn)
979               DO ji=nldit(jn),nleit(jn)
980                  png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
981               ENDDO
982            ENDDO
983         ENDDO
984         DO jj=1,ny_global
985            DO ji=1,nx_global
986               pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
987            ENDDO
988         ENDDO
989      ENDIF
990
991#endif
992
993      SELECT CASE ( cd_type )
994         CASE ( 'T' )
995            grid_loc=field_loc_center
996         CASE ( 'F' )                             
997            grid_loc=field_loc_NEcorner
998      END SELECT
999
1000      SELECT CASE ( NINT(psgn) )
1001         CASE ( -1 )
1002            field_type=field_type_vector
1003         CASE ( 1 )                             
1004            field_type=field_type_scalar
1005      END SELECT
1006
1007#if defined key_nemocice_decomp
1008      ! Ensure CICE halos are up to date
1009      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1010#else
1011!     D. Scatter pcg to CICE blocks (pc) + update halos
1012      CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
1013#endif
1014
1015      IF ( ln_ctl ) THEN
1016         WRITE(numout,*)'nemo2cice: ',varname,' min,max = ',MINVAL(pc(:,:,:)), MAXVAL(pc(:,:,:))
1017         CALL flush(numout)
1018      ENDIF
1019   END SUBROUTINE nemo2cice
1020
1021   SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
1022      !!---------------------------------------------------------------------
1023      !!                    ***  ROUTINE cice2nemo  ***
1024      !! ** Purpose :   Transfer field in CICE array to field in NEMO array.
1025#if defined key_nemocice_decomp
1026      !!             
1027      !!                NEMO and CICE PE sub domains are identical, hence
1028      !!                there is no need to gather or scatter data from
1029      !!                one PE configuration to another.
1030#else 
1031      !!                Automatically deal with scatter/gather between
1032      !!                different processors and blocks
1033      !! ** Method :    A. Gather CICE blocks (pc) into global array (pcg)
1034      !!                B. Map pcg into NEMO global array (png)
1035      !!                C. Scatter png into NEMO field (pn) for each processor
1036      !!                D. Ensure all haloes are filled in pn
1037#endif
1038      !!---------------------------------------------------------------------
1039
1040      CHARACTER(len=1), INTENT( in ) ::   &
1041          cd_type       ! nature of pn grid-point
1042          !             !   = T or F gridpoints
1043      REAL(wp), INTENT( in ) ::   &
1044          psgn          ! control of the sign change
1045          !             !   =-1 , the sign is modified following the type of b.c. used
1046          !             !   = 1 , no sign change
1047      REAL(wp), DIMENSION(jpi,jpj) :: pn
1048
1049#if defined key_nemocice_decomp
1050      INTEGER (int_kind) :: &
1051         field_type,        & ! id for type of field (scalar, vector, angle)
1052         grid_loc             ! id for location on horizontal grid
1053                              ! (center, NEcorner, Nface, Eface)
1054#else
1055      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
1056#endif
1057
1058      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
1059
1060      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
1061
1062
1063#if defined key_nemocice_decomp
1064
1065      SELECT CASE ( cd_type )
1066         CASE ( 'T' )
1067            grid_loc=field_loc_center
1068         CASE ( 'F' )                             
1069            grid_loc=field_loc_NEcorner
1070      END SELECT
1071
1072      SELECT CASE ( NINT(psgn) )
1073         CASE ( -1 )
1074            field_type=field_type_vector
1075         CASE ( 1 )                             
1076            field_type=field_type_scalar
1077      END SELECT
1078
1079      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1080
1081
1082      pn(:,:)=0.0
1083      DO jj=1,jpjm1
1084         DO ji=1,jpim1
1085            pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
1086         ENDDO
1087      ENDDO
1088
1089#else
1090
1091!      A. Gather CICE blocks (pc) into global array (pcg)
1092
1093      CALL gather_global(pcg, pc, 0, distrb_info)
1094
1095!     B. Map pcg into NEMO global array (png)
1096
1097! Need to make sure this is robust to changes in NEMO halo rows....
1098! (may be OK but not spent much time thinking about it)
1099! Note that non-existent pcg elements may be used below, but
1100! the lbclnk call on pn will replace these with sensible values
1101
1102      IF (nproc==0) THEN
1103         png(:,:,:)=0.0
1104         DO jn=1,jpnij
1105            DO jj=nldjt(jn),nlejt(jn)
1106               DO ji=nldit(jn),nleit(jn)
1107                  png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
1108               ENDDO
1109            ENDDO
1110         ENDDO
1111      ENDIF
1112
1113!     C. Scatter png into NEMO field (pn) for each processor
1114
1115      IF ( jpnij > 1) THEN
1116         CALL mppsync
1117         CALL mppscatter (png,0,pn) 
1118         CALL mppsync
1119      ELSE
1120         pn(:,:)=png(:,:,1)
1121      ENDIF
1122
1123#endif
1124
1125!     D. Ensure all haloes are filled in pn
1126
1127      CALL lbc_lnk( pn , cd_type, psgn )
1128
1129   END SUBROUTINE cice2nemo
1130
1131#else
1132   !!----------------------------------------------------------------------
1133   !!   Default option           Dummy module         NO CICE sea-ice model
1134   !!----------------------------------------------------------------------
1135CONTAINS
1136
1137   SUBROUTINE sbc_ice_cice ( kt, ksbc )     ! Dummy routine
1138      WRITE(*,*) 'sbc_ice_cice: You should not have seen this print! error?', kt
1139   END SUBROUTINE sbc_ice_cice
1140
1141   SUBROUTINE cice_sbc_init (ksbc)    ! Dummy routine
1142      WRITE(*,*) 'cice_sbc_init: You should not have seen this print! error?'
1143   END SUBROUTINE cice_sbc_init
1144
1145   SUBROUTINE cice_sbc_final     ! Dummy routine
1146      WRITE(*,*) 'cice_sbc_final: You should not have seen this print! error?'
1147   END SUBROUTINE cice_sbc_final
1148
1149#endif
1150
1151   !!======================================================================
1152END MODULE sbcice_cice
Note: See TracBrowser for help on using the repository browser.