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

source: branches/2011/dev_r2802_UKMO8_cice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 3040

Last change on this file since 3040 was 3040, checked in by charris, 13 years ago

#665 Tidied up comments, history information in headers etc.

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