source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_sit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11939

Last change on this file since 11939 was 11939, checked in by dcarneir, 10 months ago

Changing GO6 package to include sea ice thickness DA

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