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

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

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

Last change on this file since 8899 was 8899, checked in by hadlh, 6 years ago

Removed unecessary print statements and stopped the passing back of salf and freshwater fluxes as these are output in CICE now.

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