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_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 9 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

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