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 NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 10764

Last change on this file since 10764 was 10764, checked in by jcastill, 5 years ago

Minimal set of changes for coupling order

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