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

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

source: branches/UKMO/dev_3841_sbc/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 4827

Last change on this file since 4827 was 4827, checked in by charris, 9 years ago

Some demonstration code changes.

File size: 35.7 KB
Line 
1MODULE sbcice_cice
2   !!======================================================================
3   !!                       ***  MODULE  sbcice_cice  ***
4   !! To couple with sea ice model CICE (LANL)
5   !!=====================================================================
6#if defined key_cice
7   !!----------------------------------------------------------------------
8   !!   'key_cice' :                                     CICE sea-ice model
9   !!----------------------------------------------------------------------
10   !!   sbc_ice_cice  : sea-ice model time-stepping and update ocean sbc over ice-covered area
11   !!   
12   !!   
13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE domvvl
17   USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic
18   USE in_out_manager  ! I/O manager
19   USE lib_mpp         ! distributed memory computing library
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21   USE wrk_nemo        ! work arrays
22   USE timing          ! Timing
23   USE daymod          ! calendar
24   USE fld_def
25
26   USE sbc_oce         ! Surface boundary condition: ocean fields
27   USE sbc_ice         ! Surface boundary condition: ice   fields
28   USE sbcblk_core     ! Surface boundary condition: CORE bulk
29   USE sbccpl
30   USE sbcget
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   USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
41   USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
42                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm,     &
43                fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt,          &
44                flatn_f,fsurfn_f,fcondtopn_f,                    &
45                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
46                swvdr,swvdf,swidr,swidf
47   USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
48   USE ice_atmo, only: calc_strair
49   USE ice_therm_vertical, only: calc_Tsfc
50
51   USE CICE_InitMod
52   USE CICE_RunMod
53   USE CICE_FinalMod
54
55   IMPLICIT NONE
56   PRIVATE
57
58   !! * Routine accessibility
59   PUBLIC cice_sbc_init   ! routine called by sbc_init
60   PUBLIC cice_sbc_final  ! routine called by sbc_final
61   PUBLIC sbc_ice_cice    ! routine called by sbc
62
63   INTEGER , PARAMETER ::   ji_off = INT ( (jpiglo - nx_global) / 2 )
64   INTEGER , PARAMETER ::   jj_off = INT ( (jpjglo - ny_global) / 2 )
65
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE ::   png     ! local array used in sbc_cice_ice
67
68   !! * Substitutions
69#  include "domzgr_substitute.h90"
70
71CONTAINS
72
73   INTEGER FUNCTION sbc_ice_cice_alloc()
74      !!----------------------------------------------------------------------
75      !!                ***  FUNCTION sbc_ice_cice_alloc  ***
76      !!----------------------------------------------------------------------
77      ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
78      IF( lk_mpp                 )   CALL mpp_sum ( sbc_ice_cice_alloc )
79      IF( sbc_ice_cice_alloc > 0 )   CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
80   END FUNCTION sbc_ice_cice_alloc
81
82   SUBROUTINE sbc_ice_cice( kt, nsbc )
83      !!---------------------------------------------------------------------
84      !!                  ***  ROUTINE sbc_ice_cice  ***
85      !!                   
86      !! ** Purpose :   update the ocean surface boundary condition via the
87      !!                CICE Sea Ice Model time stepping
88      !!
89      !! ** Method  : - Get any extra forcing fields for CICE 
90      !!              - Prepare forcing fields
91      !!              - CICE model time stepping
92      !!              - call the routine that computes mass and
93      !!                heat fluxes at the ice/ocean interface
94      !!
95      !! ** Action  : - time evolution of the CICE sea-ice model
96      !!              - update all sbc variables below sea-ice:
97      !!                utau, vtau, qns , qsr, emp , sfx
98      !!---------------------------------------------------------------------
99      INTEGER, INTENT(in) ::   kt      ! ocean time step
100      INTEGER, INTENT(in) ::   nsbc    ! surface forcing type
101      !!----------------------------------------------------------------------
102      !
103      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_cice')
104      !
105      !                                        !----------------------!
106      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
107         !                                     !----------------------!
108
109         ! Make sure any fluxes required for CICE are set
110         IF ( nsbc == 2 )  THEN
111            CALL cice_sbc_force(kt)
112         ELSE IF ( nsbc == 5 ) THEN
113            CALL sbc_cpl_ice_flx( 1.0-fr_i  )
114         ENDIF
115
116         CALL cice_sbc_in ( kt, nsbc )
117         CALL CICE_Run
118         CALL cice_sbc_out ( kt, nsbc )
119
120         IF ( nsbc == 5 )  CALL cice_sbc_hadgam(kt+1)
121
122      ENDIF                                          ! End sea-ice time step only
123      !
124      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_cice')
125
126   END SUBROUTINE sbc_ice_cice
127
128   SUBROUTINE cice_sbc_init (nsbc)
129      !!---------------------------------------------------------------------
130      !!                    ***  ROUTINE cice_sbc_init  ***
131      !! ** Purpose: Initialise ice related fields for NEMO and coupling
132      !!
133      INTEGER, INTENT( in  ) ::   nsbc                ! surface forcing type
134      REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
135      REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar
136      INTEGER  ::   ji, jj, jl                        ! dummy loop indices
137      !!---------------------------------------------------------------------
138
139      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_init')
140      !
141      CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
142      !
143      IF(lwp) WRITE(numout,*)'cice_sbc_init'
144
145! Initialize CICE
146      CALL CICE_Initialize
147
148! Do some CICE consistency checks
149      IF ( (nsbc == 2) .OR. (nsbc == 5) ) THEN
150         IF ( calc_strair .OR. calc_Tsfc ) THEN
151            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
152         ENDIF
153      ELSEIF (nsbc == 4) THEN
154         IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
155            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
156         ENDIF
157      ENDIF
158
159
160! allocate sbc_ice and sbc_cice arrays
161      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' )
162      IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
163
164! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
165      IF( .NOT. ln_rstart ) THEN
166         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
167         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
168      ENDIF
169
170      fr_iu(:,:)=0.0
171      fr_iv(:,:)=0.0
172
173      CALL cice2nemo(aice,fr_i, 'T', 1. )
174      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
175         DO jl=1,ncat
176            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
177         ENDDO
178      ENDIF
179
180! T point to U point
181! T point to V point
182      DO jj=1,jpjm1
183         DO ji=1,jpim1
184            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
185            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
186         ENDDO
187      ENDDO
188
189      CALL lbc_lnk ( fr_iu , 'U', 1. )
190      CALL lbc_lnk ( fr_iv , 'V', 1. )
191
192      !                                      ! embedded sea ice
193      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
194         CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
195         CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
196         snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
197         snwice_mass_b(:,:) = snwice_mass(:,:)
198      ELSE
199         snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges
200         snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
201      ENDIF
202      IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart :
203         &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh belew sea-ice area
204         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
205         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
206         !
207         ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n]
208         !       which were previously set in domvvl
209         IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only???
210            DO jj = 1, jpjm1
211               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
212                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
213                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
214                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
215                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
216                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
217                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
218                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
219                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
220                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
221                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
222                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
223               END DO
224            END DO
225            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
226            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
227            DO jj = 1, jpjm1
228               DO ji = 1, jpim1      ! NO Vector Opt.
229                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
230                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
231                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
232                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
233               END DO
234            END DO
235            CALL lbc_lnk( sshf_n, 'F', 1. )
236          ENDIF
237      ENDIF
238 
239      CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
240      !
241      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init')
242      !
243   END SUBROUTINE cice_sbc_init
244
245   
246   SUBROUTINE cice_sbc_in (kt, nsbc)
247      !!---------------------------------------------------------------------
248      !!                    ***  ROUTINE cice_sbc_in  ***
249      !! ** Purpose: Set coupling fields and pass to CICE
250      !!---------------------------------------------------------------------
251      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
252      INTEGER, INTENT(in   ) ::   nsbc ! surface forcing type
253
254      INTEGER  ::   ji, jj, jl                   ! dummy loop indices     
255      REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice
256      REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
257      REAL(wp) ::   zintb, zintn  ! dummy argument
258      !!---------------------------------------------------------------------
259
260      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_in')
261      !
262      CALL wrk_alloc( jpi,jpj, ztmp, zpice )
263      CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
264
265      IF( kt == nit000 )  THEN
266         IF(lwp) WRITE(numout,*)'cice_sbc_in'
267      ENDIF
268
269      ztmp(:,:)=0.0
270
271! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
272! the first time-step)
273
274! forced and coupled case
275
276      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
277
278         ztmpn(:,:,:)=0.0
279
280! x comp of wind stress (CI_1)
281! U point to F point
282         DO jj=1,jpjm1
283            DO ji=1,jpi
284               ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
285                                    + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
286            ENDDO
287         ENDDO
288         CALL nemo2cice(ztmp,strax,'F', -1. )
289
290! y comp of wind stress (CI_2)
291! V point to F point
292         DO jj=1,jpj
293            DO ji=1,jpim1
294               ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
295                                    + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
296            ENDDO
297         ENDDO
298         CALL nemo2cice(ztmp,stray,'F', -1. )
299
300! Surface downward latent heat flux (CI_5)
301         IF (nsbc == 2) THEN
302            DO jl=1,ncat
303               ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
304            ENDDO
305         ELSE
306! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
307            qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub
308! End of temporary code
309            DO jj=1,jpj
310               DO ji=1,jpi
311                  IF (fr_i(ji,jj).eq.0.0) THEN
312                     DO jl=1,ncat
313                        ztmpn(ji,jj,jl)=0.0
314                     ENDDO
315                     ! This will then be conserved in CICE
316                     ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
317                  ELSE
318                     DO jl=1,ncat
319                        ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
320                     ENDDO
321                  ENDIF
322               ENDDO
323            ENDDO
324         ENDIF
325         DO jl=1,ncat
326            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
327
328! GBM conductive flux through ice (CI_6)
329!  Convert to GBM
330            IF (nsbc == 2) THEN
331               ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
332            ELSE
333               ztmp(:,:) = botmelt(:,:,jl)
334            ENDIF
335            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
336
337! GBM surface heat flux (CI_7)
338!  Convert to GBM
339            IF (nsbc == 2) THEN
340               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 
341            ELSE
342               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
343            ENDIF
344            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
345         ENDDO
346
347      ELSE IF (nsbc == 4) THEN
348
349! Pass CORE forcing fields to CICE (which will calculate heat fluxes etc itself)
350! x comp and y comp of atmosphere surface wind (CICE expects on T points)
351         ztmp(:,:) = wndi_ice(:,:)
352         CALL nemo2cice(ztmp,uatm,'T', -1. )
353         ztmp(:,:) = wndj_ice(:,:)
354         CALL nemo2cice(ztmp,vatm,'T', -1. )
355         ztmp(:,:) = SQRT ( wndi_ice(:,:)**2 + wndj_ice(:,:)**2 )
356         CALL nemo2cice(ztmp,wind,'T', 1. )    ! Wind speed (m/s)
357         ztmp(:,:) = qsr_ice(:,:,1)
358         CALL nemo2cice(ztmp,fsw,'T', 1. )     ! Incoming short-wave (W/m^2)
359         ztmp(:,:) = qlw_ice(:,:,1)
360         CALL nemo2cice(ztmp,flw,'T', 1. )     ! Incoming long-wave (W/m^2)
361         ztmp(:,:) = tatm_ice(:,:)
362         CALL nemo2cice(ztmp,Tair,'T', 1. )    ! Air temperature (K)
363         CALL nemo2cice(ztmp,potT,'T', 1. )    ! Potential temp (K)
364! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows 
365         ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )   
366                                               ! Constant (101000.) atm pressure assumed
367         CALL nemo2cice(ztmp,rhoa,'T', 1. )    ! Air density (kg/m^3)
368         ztmp(:,:) = qatm_ice(:,:)
369         CALL nemo2cice(ztmp,Qa,'T', 1. )      ! Specific humidity (kg/kg)
370         ztmp(:,:)=10.0
371         CALL nemo2cice(ztmp,zlvl,'T', 1. )    ! Atmos level height (m)
372
373! May want to check all values are physically realistic (as in CICE routine
374! prepare_forcing)?
375
376! Divide shortwave into spectral bands (as in prepare_forcing)
377         ztmp(:,:)=qsr_ice(:,:,1)*frcvdr       ! visible direct
378         CALL nemo2cice(ztmp,swvdr,'T', 1. )             
379         ztmp(:,:)=qsr_ice(:,:,1)*frcvdf       ! visible diffuse
380         CALL nemo2cice(ztmp,swvdf,'T', 1. )             
381         ztmp(:,:)=qsr_ice(:,:,1)*frcidr       ! near IR direct
382         CALL nemo2cice(ztmp,swidr,'T', 1. )
383         ztmp(:,:)=qsr_ice(:,:,1)*frcidf       ! near IR diffuse
384         CALL nemo2cice(ztmp,swidf,'T', 1. )
385
386      ENDIF
387
388! Snowfall
389! Ensure fsnow is positive (as in CICE routine prepare_forcing) 
390      ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0) 
391      CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
392
393! Rainfall
394      ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
395      CALL nemo2cice(ztmp,frain,'T', 1. ) 
396
397! Freezing/melting potential
398! Calculated over NEMO leapfrog timestep (hence 2*dt)
399      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:,1)*(Tocnfrz-sst_m(:,:))/(2.0*dt)
400
401      ztmp(:,:) = nfrzmlt(:,:)
402      CALL nemo2cice(ztmp,frzmlt,'T', 1. )
403
404! SST  and SSS
405
406      CALL nemo2cice(sst_m,sst,'T', 1. )
407      CALL nemo2cice(sss_m,sss,'T', 1. )
408
409! x comp and y comp of surface ocean current
410! U point to F point
411      DO jj=1,jpjm1
412         DO ji=1,jpi
413            ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
414         ENDDO
415      ENDDO
416      CALL nemo2cice(ztmp,uocn,'F', -1. )
417
418! V point to F point
419      DO jj=1,jpj
420         DO ji=1,jpim1
421            ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
422         ENDDO
423      ENDDO
424      CALL nemo2cice(ztmp,vocn,'F', -1. )
425
426      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
427          !
428          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
429          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
430         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
431          !
432          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
433          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
434         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
435          !
436         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
437          !
438         !
439      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
440         zpice(:,:) = ssh_m(:,:)
441      ENDIF
442
443! x comp and y comp of sea surface slope (on F points)
444! T point to F point
445      DO jj=1,jpjm1
446         DO ji=1,jpim1
447            ztmp(ji,jj)=0.5 * (  (zpice(ji+1,jj  )-zpice(ji,jj  ))/e1u(ji,jj  )   &
448                               + (zpice(ji+1,jj+1)-zpice(ji,jj+1))/e1u(ji,jj+1) ) & 
449                            *  fmask(ji,jj,1)
450         ENDDO
451      ENDDO
452      CALL nemo2cice(ztmp,ss_tltx,'F', -1. )
453
454! T point to F point
455      DO jj=1,jpjm1
456         DO ji=1,jpim1
457            ztmp(ji,jj)=0.5 * (  (zpice(ji  ,jj+1)-zpice(ji  ,jj))/e2v(ji  ,jj)   &
458                               + (zpice(ji+1,jj+1)-zpice(ji+1,jj))/e2v(ji+1,jj) ) &
459                            *  fmask(ji,jj,1)
460         ENDDO
461      ENDDO
462      CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
463
464      CALL wrk_dealloc( jpi,jpj, ztmp )
465      CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
466      !
467      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_in')
468      !
469   END SUBROUTINE cice_sbc_in
470
471
472   SUBROUTINE cice_sbc_out (kt,nsbc)
473      !!---------------------------------------------------------------------
474      !!                    ***  ROUTINE cice_sbc_out  ***
475      !! ** Purpose: Get fields from CICE and set surface fields for NEMO
476      !!---------------------------------------------------------------------
477      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
478      INTEGER, INTENT( in  ) ::   nsbc ! surface forcing type
479     
480      INTEGER  ::   ji, jj, jl                 ! dummy loop indices
481      REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
482      !!---------------------------------------------------------------------
483
484      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_out')
485      !
486      CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
487     
488      IF( kt == nit000 )  THEN
489         IF(lwp) WRITE(numout,*)'cice_sbc_out'
490      ENDIF
491     
492! x comp of ocean-ice stress
493      CALL cice2nemo(strocnx,ztmp1,'F', -1. )
494      ss_iou(:,:)=0.0
495! F point to U point
496      DO jj=2,jpjm1
497         DO ji=2,jpim1
498            ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
499         ENDDO
500      ENDDO
501      CALL lbc_lnk( ss_iou , 'U', -1. )
502
503! y comp of ocean-ice stress
504      CALL cice2nemo(strocny,ztmp1,'F', -1. )
505      ss_iov(:,:)=0.0
506! F point to V point
507
508      DO jj=1,jpjm1
509         DO ji=2,jpim1
510            ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
511         ENDDO
512      ENDDO
513      CALL lbc_lnk( ss_iov , 'V', -1. )
514
515! x and y comps of surface stress
516! Combine wind stress and ocean-ice stress
517! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
518
519      utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
520      vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)     
521
522! Freshwater fluxes
523
524      IF (nsbc == 2) THEN
525! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
526! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below
527! Not ideal since aice won't be the same as in the atmosphere. 
528! Better to use evap and tprecip? (but for now don't read in evap in this case)
529         emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
530      ELSE IF (nsbc == 4) THEN
531         emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)       
532      ELSE IF (nsbc ==5) THEN
533! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
534! This is currently as required with the coupling fields from the UM atmosphere
535         emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:) 
536      ENDIF
537
538      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
539      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
540
541! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
542! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
543! which has to be compensated for by the ocean salinity potentially going negative
544! This check breaks conservation but seems reasonable until we have prognostic ice salinity
545! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
546      WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
547      sfx(:,:)=ztmp2(:,:)*1000.0
548      emp(:,:)=emp(:,:)-ztmp1(:,:)
549 
550      CALL lbc_lnk( emp , 'T', 1. )
551      CALL lbc_lnk( sfx , 'T', 1. )
552
553! Solar penetrative radiation and non solar surface heat flux
554
555! Scale qsr and qns according to ice fraction (bulk formulae only)
556
557      IF (nsbc == 4) THEN
558         qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
559         qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
560      ENDIF
561! Take into account snow melting except for fully coupled when already in qns_tot
562      IF (nsbc == 5) THEN
563         qsr(:,:)= qsr_tot(:,:)
564         qns(:,:)= qns_tot(:,:)
565      ELSE
566         qns(:,:)= qns(:,:)-sprecip(:,:)*Lfresh*(1.0-fr_i(:,:))
567      ENDIF
568
569! Now add in ice / snow related terms
570! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
571      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
572      qsr(:,:)=qsr(:,:)+ztmp1(:,:)
573      CALL lbc_lnk( qsr , 'T', 1. )
574
575      DO jj=1,jpj
576         DO ji=1,jpi
577            nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
578         ENDDO
579      ENDDO
580
581      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
582      qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
583
584      CALL lbc_lnk( qns , 'T', 1. )
585
586! Prepare for the following CICE time-step
587
588      CALL cice2nemo(aice,fr_i,'T', 1. )
589      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
590         DO jl=1,ncat
591            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
592         ENDDO
593      ENDIF
594
595! T point to U point
596! T point to V point
597      DO jj=1,jpjm1
598         DO ji=1,jpim1
599            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
600            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
601         ENDDO
602      ENDDO
603
604      CALL lbc_lnk ( fr_iu , 'U', 1. )
605      CALL lbc_lnk ( fr_iv , 'V', 1. )
606
607      !                                      ! embedded sea ice
608      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
609         CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
610         CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
611         snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
612         snwice_mass_b(:,:) = snwice_mass(:,:)
613         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
614      ENDIF
615
616! Release work space
617
618      CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
619      !
620      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_out')
621      !
622   END SUBROUTINE cice_sbc_out
623
624
625#if defined key_oasis3 || defined key_oasis4
626   SUBROUTINE cice_sbc_hadgam( kt )
627      !!---------------------------------------------------------------------
628      !!                    ***  ROUTINE cice_sbc_hadgam  ***
629      !! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
630      !!
631      !!
632      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
633      !!---------------------------------------------------------------------
634
635      INTEGER  ::   jl                        ! dummy loop index
636      INTEGER  ::   ierror
637
638      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_hadgam')
639      !
640      IF( kt == nit000 )  THEN
641         IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
642!         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
643      ENDIF
644
645      !                                         ! =========================== !
646      !                                         !   Prepare Coupling fields   !
647      !                                         ! =========================== !
648
649! x and y comp of ice velocity
650
651      CALL cice2nemo(uvel,u_ice,'F', -1. )
652      CALL cice2nemo(vvel,v_ice,'F', -1. )
653
654! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out 
655
656! Snow and ice thicknesses (CO_2 and CO_3)
657
658      DO jl = 1,ncat
659         CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. )
660         CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
661      ENDDO
662      !
663      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_hadgam')
664      !
665   END SUBROUTINE cice_sbc_hadgam
666
667#else
668   SUBROUTINE cice_sbc_hadgam( kt )    ! Dummy routine
669      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
670      WRITE(*,*) 'cice_sbc_hadgam: You should not have seen this print! error?'
671   END SUBROUTINE cice_sbc_hadgam
672#endif
673
674   SUBROUTINE cice_sbc_final
675      !!---------------------------------------------------------------------
676      !!                    ***  ROUTINE cice_sbc_final  ***
677      !! ** Purpose: Finalize CICE
678      !!---------------------------------------------------------------------
679
680      IF(lwp) WRITE(numout,*)'cice_sbc_final'
681
682      CALL CICE_Finalize
683
684   END SUBROUTINE cice_sbc_final
685
686   SUBROUTINE cice_sbc_force (kt)
687      !!---------------------------------------------------------------------
688      !!                    ***  ROUTINE cice_sbc_force  ***
689      !! ** Purpose : Provide CICE forcing from files
690      !!
691      !!---------------------------------------------------------------------
692      !! ** Method  :   Set forcing fields
693      !!     
694      !!  snowfall   
695      !!  rainfall   
696      !!  sublimation rate   
697      !!  topmelt (category)
698      !!  botmelt (category)
699      !!
700      !! History :
701      !!----------------------------------------------------------------------
702      !! * Modules used
703
704      INTEGER, INTENT( in  ) ::   kt ! ocean time step
705
706      ! Assume the fluxes have already been obtained somewhere....
707
708      ! set the fluxes from read fields
709      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
710      tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
711! May be better to do this conversion somewhere else
712      qla_ice(:,:,1) = -Lsub*sf(jp_ievp)%fnow(:,:,1)
713      topmelt(:,:,:) = sf(jp_topm)%fnow(:,:,:)
714      botmelt(:,:,:) = sf(jp_botm)%fnow(:,:,:)
715
716   END SUBROUTINE cice_sbc_force
717
718   SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
719      !!---------------------------------------------------------------------
720      !!                    ***  ROUTINE nemo2cice  ***
721      !! ** Purpose :   Transfer field in NEMO array to field in CICE array. 
722#if defined key_nemocice_decomp
723      !!             
724      !!                NEMO and CICE PE sub domains are identical, hence
725      !!                there is no need to gather or scatter data from
726      !!                one PE configuration to another.
727#else
728      !!                Automatically gather/scatter between
729      !!                different processors and blocks
730      !! ** Method :    A. Ensure all haloes are filled in NEMO field (pn)
731      !!                B. Gather pn into global array (png)
732      !!                C. Map png into CICE global array (pcg)
733      !!                D. Scatter pcg to CICE blocks (pc) + update haloes 
734#endif
735      !!---------------------------------------------------------------------
736
737      CHARACTER(len=1), INTENT( in ) ::   &
738          cd_type       ! nature of pn grid-point
739          !             !   = T or F gridpoints
740      REAL(wp), INTENT( in ) ::   &
741          psgn          ! control of the sign change
742          !             !   =-1 , the sign is modified following the type of b.c. used
743          !             !   = 1 , no sign change
744      REAL(wp), DIMENSION(jpi,jpj) :: pn
745#if !defined key_nemocice_decomp
746      REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
747      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
748#endif
749      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
750      INTEGER (int_kind) :: &
751         field_type,        &! id for type of field (scalar, vector, angle)
752         grid_loc            ! id for location on horizontal grid
753                            !  (center, NEcorner, Nface, Eface)
754
755      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
756
757!     A. Ensure all haloes are filled in NEMO field (pn)
758
759      CALL lbc_lnk( pn , cd_type, psgn )
760
761#if defined key_nemocice_decomp
762
763      ! Copy local domain data from NEMO to CICE field
764      pc(:,:,1)=0.0
765      DO jj=2,ny_block-1
766         DO ji=2,nx_block-1
767            pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
768         ENDDO
769      ENDDO
770
771#else
772
773!     B. Gather pn into global array (png)
774
775      IF ( jpnij > 1) THEN
776         CALL mppsync
777         CALL mppgather (pn,0,png) 
778         CALL mppsync
779      ELSE
780         png(:,:,1)=pn(:,:)
781      ENDIF
782
783!     C. Map png into CICE global array (pcg)
784
785! Need to make sure this is robust to changes in NEMO halo rows....
786! (may be OK but not 100% sure)
787
788      IF (nproc==0) THEN     
789!        pcg(:,:)=0.0
790         DO jn=1,jpnij
791            DO jj=nldjt(jn),nlejt(jn)
792               DO ji=nldit(jn),nleit(jn)
793                  png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
794               ENDDO
795            ENDDO
796         ENDDO
797         DO jj=1,ny_global
798            DO ji=1,nx_global
799               pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
800            ENDDO
801         ENDDO
802      ENDIF
803
804#endif
805
806      SELECT CASE ( cd_type )
807         CASE ( 'T' )
808            grid_loc=field_loc_center
809         CASE ( 'F' )                             
810            grid_loc=field_loc_NEcorner
811      END SELECT
812
813      SELECT CASE ( NINT(psgn) )
814         CASE ( -1 )
815            field_type=field_type_vector
816         CASE ( 1 )                             
817            field_type=field_type_scalar
818      END SELECT
819
820#if defined key_nemocice_decomp
821      ! Ensure CICE halos are up to date
822      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
823#else
824!     D. Scatter pcg to CICE blocks (pc) + update halos
825      CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
826#endif
827
828   END SUBROUTINE nemo2cice
829
830   SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
831      !!---------------------------------------------------------------------
832      !!                    ***  ROUTINE cice2nemo  ***
833      !! ** Purpose :   Transfer field in CICE array to field in NEMO array.
834#if defined key_nemocice_decomp
835      !!             
836      !!                NEMO and CICE PE sub domains are identical, hence
837      !!                there is no need to gather or scatter data from
838      !!                one PE configuration to another.
839#else 
840      !!                Automatically deal with scatter/gather between
841      !!                different processors and blocks
842      !! ** Method :    A. Gather CICE blocks (pc) into global array (pcg)
843      !!                B. Map pcg into NEMO global array (png)
844      !!                C. Scatter png into NEMO field (pn) for each processor
845      !!                D. Ensure all haloes are filled in pn
846#endif
847      !!---------------------------------------------------------------------
848
849      CHARACTER(len=1), INTENT( in ) ::   &
850          cd_type       ! nature of pn grid-point
851          !             !   = T or F gridpoints
852      REAL(wp), INTENT( in ) ::   &
853          psgn          ! control of the sign change
854          !             !   =-1 , the sign is modified following the type of b.c. used
855          !             !   = 1 , no sign change
856      REAL(wp), DIMENSION(jpi,jpj) :: pn
857
858#if defined key_nemocice_decomp
859      INTEGER (int_kind) :: &
860         field_type,        & ! id for type of field (scalar, vector, angle)
861         grid_loc             ! id for location on horizontal grid
862                              ! (center, NEcorner, Nface, Eface)
863#else
864      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
865#endif
866
867      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
868
869      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
870
871
872#if defined key_nemocice_decomp
873
874      SELECT CASE ( cd_type )
875         CASE ( 'T' )
876            grid_loc=field_loc_center
877         CASE ( 'F' )                             
878            grid_loc=field_loc_NEcorner
879      END SELECT
880
881      SELECT CASE ( NINT(psgn) )
882         CASE ( -1 )
883            field_type=field_type_vector
884         CASE ( 1 )                             
885            field_type=field_type_scalar
886      END SELECT
887
888      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
889
890
891      pn(:,:)=0.0
892      DO jj=1,jpjm1
893         DO ji=1,jpim1
894            pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
895         ENDDO
896      ENDDO
897
898#else
899
900!      A. Gather CICE blocks (pc) into global array (pcg)
901
902      CALL gather_global(pcg, pc, 0, distrb_info)
903
904!     B. Map pcg into NEMO global array (png)
905
906! Need to make sure this is robust to changes in NEMO halo rows....
907! (may be OK but not spent much time thinking about it)
908! Note that non-existent pcg elements may be used below, but
909! the lbclnk call on pn will replace these with sensible values
910
911      IF (nproc==0) THEN
912         png(:,:,:)=0.0
913         DO jn=1,jpnij
914            DO jj=nldjt(jn),nlejt(jn)
915               DO ji=nldit(jn),nleit(jn)
916                  png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
917               ENDDO
918            ENDDO
919         ENDDO
920      ENDIF
921
922!     C. Scatter png into NEMO field (pn) for each processor
923
924      IF ( jpnij > 1) THEN
925         CALL mppsync
926         CALL mppscatter (png,0,pn) 
927         CALL mppsync
928      ELSE
929         pn(:,:)=png(:,:,1)
930      ENDIF
931
932#endif
933
934!     D. Ensure all haloes are filled in pn
935
936      CALL lbc_lnk( pn , cd_type, psgn )
937
938   END SUBROUTINE cice2nemo
939
940#else
941   !!----------------------------------------------------------------------
942   !!   Default option           Dummy module         NO CICE sea-ice model
943   !!----------------------------------------------------------------------
944CONTAINS
945
946   SUBROUTINE sbc_ice_cice ( kt, nsbc )     ! Dummy routine
947      WRITE(*,*) 'sbc_ice_cice: You should not have seen this print! error?', kt
948   END SUBROUTINE sbc_ice_cice
949
950   SUBROUTINE cice_sbc_init (nsbc)    ! Dummy routine
951      WRITE(*,*) 'cice_sbc_init: You should not have seen this print! error?'
952   END SUBROUTINE cice_sbc_init
953
954   SUBROUTINE cice_sbc_final     ! Dummy routine
955      WRITE(*,*) 'cice_sbc_final: You should not have seen this print! error?'
956   END SUBROUTINE cice_sbc_final
957
958#endif
959
960   !!======================================================================
961END MODULE sbcice_cice
Note: See TracBrowser for help on using the repository browser.