source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 10394

Last change on this file since 10394 was 10394, checked in by jcastill, 23 months ago

Merged r6232_hadgem3_cplseq@7460 branch

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