source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5516

Last change on this file since 5516 was 5516, checked in by smasson, 5 years ago

work_alloc/work_dealloc errors

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