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

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

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

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

Bug fix

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