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/GO6_package_r8356_CICEnudging/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/GO6_package_r8356_CICEnudging/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8618

Last change on this file since 8618 was 8618, checked in by hadlh, 7 years ago

Changes for DePreSys3 (coupled) sea-ice nudging, including code to read in ancillary file and pass to CICE.

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