source: branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8190

Last change on this file since 8190 was 8190, checked in by jwhile, 3 years ago

Update due to Tim's review

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