source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8400

Last change on this file since 8400 was 8400, checked in by timgraham, 3 years ago

GMED ticket 335:

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