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

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

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5030

Last change on this file since 5030 was 5030, checked in by timgraham, 9 years ago

Added changes for coupling of meltpond fraction and depth between NEMO-CICE and atmosphere model

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