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

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

source: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 6507

Last change on this file since 6507 was 6507, checked in by timgraham, 8 years ago

First attempt at merging in science changes from GO6 package branch at v3.6 stable (Note-namelists not yet dealt with)

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