New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcice_cice.F90 in branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5591

Last change on this file since 5591 was 5591, checked in by davestorkey, 9 years ago

Bug fixes for UKMO/dev_r5107_hadgem3_cplseq branch.

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