source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11107

Last change on this file since 11107 was 11107, checked in by frrh, 18 months ago

Commit changes from Dan Copsey's sea ice heat coupling
flux fixes in branch:
branches/UKMO/dev_r5518_GO6_fix_zemp_ice_10681
revisions 11028:11088.

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