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

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11935

Last change on this file since 11935 was 11935, checked in by dcarneir, 4 years ago

Changing SBC scripts (sbc_oce and sbcice_cice) to include sea ice thickness

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