source: NEMO/branches/2020/ticket2487/src/OCE/SBC/sbcice_cice.F90

Last change on this file was 13299, checked in by smueller, 3 months ago

Addition of a namelist parameter to enable the selectability of sea-level compensation for non-zero sea-ice/snow mass (see ticket #2487)

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