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 @ 5445

Last change on this file since 5445 was 5445, checked in by davestorkey, 9 years ago

Clear SVN keywords from 2015/dev_r5021_UKMO1_CICE_coupling branch.

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