Ticket #975: sbccpl.F90

File sbccpl.F90, 95.4 KB (added by kmogensen, 9 years ago)
Line 
1MODULE sbccpl
2   !!======================================================================
3   !!                       ***  MODULE  sbccpl  ***
4   !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode
5   !!======================================================================
6   !! History :  2.0  ! 2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
7   !!            3.0  ! 2008-02  (G. Madec, C Talandier)  surface module
8   !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
9   !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
10   !!----------------------------------------------------------------------
11#if defined key_oasis3 || defined key_oasis4
12   !!----------------------------------------------------------------------
13   !!   'key_oasis3' or 'key_oasis4'   Coupled Ocean/Atmosphere formulation
14   !!----------------------------------------------------------------------
15   !!   namsbc_cpl      : coupled formulation namlist
16   !!   sbc_cpl_init    : initialisation of the coupled exchanges
17   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
18   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
19   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
20   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
21   !!   sbc_cpl_snd     : send     fields to the atmosphere
22   !!----------------------------------------------------------------------
23   USE dom_oce         ! ocean space and time domain
24   USE sbc_oce         ! Surface boundary condition: ocean fields
25   USE sbc_ice         ! Surface boundary condition: ice fields
26   USE sbcdcy          ! surface boundary condition: diurnal cycle
27   USE phycst          ! physical constants
28#if defined key_lim3
29   USE par_ice         ! ice parameters
30   USE ice             ! ice variables
31#endif
32#if defined key_lim2
33   USE par_ice_2       ! ice parameters
34   USE ice_2           ! ice variables
35#endif
36#if defined key_oasis3
37   USE cpl_oasis3      ! OASIS3 coupling
38#endif
39#if defined key_oasis4
40   USE cpl_oasis4      ! OASIS4 coupling
41#endif
42   USE geo2ocean       !
43   USE restart         !
44   USE oce   , ONLY : tsn, un, vn
45   USE albedo          !
46   USE in_out_manager  ! I/O manager
47   USE iom             ! NetCDF library
48   USE lib_mpp         ! distribued memory computing library
49   USE wrk_nemo        ! work arrays
50   USE timing          ! Timing
51   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
52#if defined key_cpl_carbon_cycle
53   USE p4zflx, ONLY : oce_co2
54#endif
55   USE diaar5, ONLY :   lk_diaar5
56#if defined key_cice
57   USE ice_domain_size, only: ncat
58#endif
59   IMPLICIT NONE
60   PRIVATE
61
62   PUBLIC   sbc_cpl_rcv        ! routine called by sbc_ice_lim(_2).F90
63   PUBLIC   sbc_cpl_snd        ! routine called by step.F90
64   PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
65   PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).F90
66
67   INTEGER, PARAMETER ::   jpr_otx1   =  1            ! 3 atmosphere-ocean stress components on grid 1
68   INTEGER, PARAMETER ::   jpr_oty1   =  2            !
69   INTEGER, PARAMETER ::   jpr_otz1   =  3            !
70   INTEGER, PARAMETER ::   jpr_otx2   =  4            ! 3 atmosphere-ocean stress components on grid 2
71   INTEGER, PARAMETER ::   jpr_oty2   =  5            !
72   INTEGER, PARAMETER ::   jpr_otz2   =  6            !
73   INTEGER, PARAMETER ::   jpr_itx1   =  7            ! 3 atmosphere-ice   stress components on grid 1
74   INTEGER, PARAMETER ::   jpr_ity1   =  8            !
75   INTEGER, PARAMETER ::   jpr_itz1   =  9            !
76   INTEGER, PARAMETER ::   jpr_itx2   = 10            ! 3 atmosphere-ice   stress components on grid 2
77   INTEGER, PARAMETER ::   jpr_ity2   = 11            !
78   INTEGER, PARAMETER ::   jpr_itz2   = 12            !
79   INTEGER, PARAMETER ::   jpr_qsroce = 13            ! Qsr above the ocean
80   INTEGER, PARAMETER ::   jpr_qsrice = 14            ! Qsr above the ice
81   INTEGER, PARAMETER ::   jpr_qsrmix = 15 
82   INTEGER, PARAMETER ::   jpr_qnsoce = 16            ! Qns above the ocean
83   INTEGER, PARAMETER ::   jpr_qnsice = 17            ! Qns above the ice
84   INTEGER, PARAMETER ::   jpr_qnsmix = 18
85   INTEGER, PARAMETER ::   jpr_rain   = 19            ! total liquid precipitation (rain)
86   INTEGER, PARAMETER ::   jpr_snow   = 20            ! solid precipitation over the ocean (snow)
87   INTEGER, PARAMETER ::   jpr_tevp   = 21            ! total evaporation
88   INTEGER, PARAMETER ::   jpr_ievp   = 22            ! solid evaporation (sublimation)
89   INTEGER, PARAMETER ::   jpr_sbpr   = 23            ! sublimation - liquid precipitation - solid precipitation
90   INTEGER, PARAMETER ::   jpr_semp   = 24            ! solid freshwater budget (sublimation - snow)
91   INTEGER, PARAMETER ::   jpr_oemp   = 25            ! ocean freshwater budget (evap - precip)
92   INTEGER, PARAMETER ::   jpr_w10m   = 26            ! 10m wind
93   INTEGER, PARAMETER ::   jpr_dqnsdt = 27            ! d(Q non solar)/d(temperature)
94   INTEGER, PARAMETER ::   jpr_rnf    = 28            ! runoffs
95   INTEGER, PARAMETER ::   jpr_cal    = 29            ! calving
96   INTEGER, PARAMETER ::   jpr_taum   = 30            ! wind stress module
97   INTEGER, PARAMETER ::   jpr_co2    = 31
98   INTEGER, PARAMETER ::   jpr_topm   = 32            ! topmeltn
99   INTEGER, PARAMETER ::   jpr_botm   = 33            ! botmeltn
100   INTEGER, PARAMETER ::   jprcv      = 33            ! total number of fields received
101
102   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction
103   INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
104   INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
105   INTEGER, PARAMETER ::   jps_tmix   =  4            ! mixed temperature (ocean+ice)
106   INTEGER, PARAMETER ::   jps_albice =  5            ! ice   albedo
107   INTEGER, PARAMETER ::   jps_albmix =  6            ! mixed albedo
108   INTEGER, PARAMETER ::   jps_hice   =  7            ! ice  thickness
109   INTEGER, PARAMETER ::   jps_hsnw   =  8            ! snow thickness
110   INTEGER, PARAMETER ::   jps_ocx1   =  9            ! ocean current on grid 1
111   INTEGER, PARAMETER ::   jps_ocy1   = 10            !
112   INTEGER, PARAMETER ::   jps_ocz1   = 11            !
113   INTEGER, PARAMETER ::   jps_ivx1   = 12            ! ice   current on grid 1
114   INTEGER, PARAMETER ::   jps_ivy1   = 13            !
115   INTEGER, PARAMETER ::   jps_ivz1   = 14            !
116   INTEGER, PARAMETER ::   jps_co2    = 15
117   INTEGER, PARAMETER ::   jpsnd      = 15            ! total number of fields sended
118
119   !                                                         !!** namelist namsbc_cpl **
120   TYPE ::   FLD_C
121      CHARACTER(len = 32) ::   cldes                  ! desciption of the coupling strategy
122      CHARACTER(len = 32) ::   clcat                  ! multiple ice categories strategy
123      CHARACTER(len = 32) ::   clvref                 ! reference of vector ('spherical' or 'cartesian')
124      CHARACTER(len = 32) ::   clvor                  ! orientation of vector fields ('eastward-northward' or 'local grid')
125      CHARACTER(len = 32) ::   clvgrd                 ! grids on which is located the vector fields
126   END TYPE FLD_C
127   ! Send to the atmosphere                           !
128   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                       
129   ! Received from the atmosphere                     !
130   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
131   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                       
132
133   TYPE ::   DYNARR     
134      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
135   END TYPE DYNARR
136
137   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                      ! all fields recieved from the atmosphere
138
139   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
140
141   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
142
143#if ! defined key_lim2   &&   ! defined key_lim3
144   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice,fr1_i0,fr2_i0          ! jpi, jpj
145   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl)
146#endif
147
148#if defined key_cice
149   INTEGER, PARAMETER ::   jpl = ncat
150#elif ! defined key_lim2   &&   ! defined key_lim3
151   INTEGER, PARAMETER ::   jpl = 1 
152   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice
153   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice
154#endif
155
156#if ! defined key_lim3   &&  ! defined key_cice
157   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i
158#endif
159
160#if ! defined key_lim3
161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  ht_i, ht_s
162#endif
163
164#if ! defined key_cice
165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  topmelt, botmelt
166#endif
167
168   !! Substitution
169#  include "vectopt_loop_substitute.h90"
170   !!----------------------------------------------------------------------
171   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
172   !! $Id: sbccpl.F90 3294 2012-01-28 16:44:18Z rblod $
173   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
174   !!----------------------------------------------------------------------
175
176CONTAINS
177 
178   INTEGER FUNCTION sbc_cpl_alloc()
179      !!----------------------------------------------------------------------
180      !!             ***  FUNCTION sbc_cpl_alloc  ***
181      !!----------------------------------------------------------------------
182      INTEGER :: ierr(4),jn
183      !!----------------------------------------------------------------------
184      ierr(:) = 0
185      !
186      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
187      !
188#if ! defined key_lim2 && ! defined key_lim3
189      ! quick patch to be able to run the coupled model without sea-ice...
190      ALLOCATE( u_ice(jpi,jpj) , fr1_i0(jpi,jpj) , tn_ice (jpi,jpj,1) ,     &
191                v_ice(jpi,jpj) , fr2_i0(jpi,jpj) , alb_ice(jpi,jpj,1),      &
192                emp_ice(jpi,jpj) , qns_ice(jpi,jpj,1) , dqns_ice(jpi,jpj,1) , STAT=ierr(2) )
193#endif
194
195#if ! defined key_lim3 && ! defined key_cice
196      ALLOCATE( a_i(jpi,jpj,jpl) , STAT=ierr(3) )
197#endif
198
199#if defined key_cice || defined key_lim2
200      ALLOCATE( ht_i(jpi,jpj,jpl) , ht_s(jpi,jpj,jpl) , STAT=ierr(4) )
201#endif
202      sbc_cpl_alloc = MAXVAL( ierr )
203      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
204      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
205      !
206   END FUNCTION sbc_cpl_alloc
207
208
209   SUBROUTINE sbc_cpl_init( k_ice )     
210      !!----------------------------------------------------------------------
211      !!             ***  ROUTINE sbc_cpl_init  ***
212      !!
213      !! ** Purpose :   Initialisation of send and recieved information from
214      !!                the atmospheric component
215      !!
216      !! ** Method  : * Read namsbc_cpl namelist
217      !!              * define the receive interface
218      !!              * define the send    interface
219      !!              * initialise the OASIS coupler
220      !!----------------------------------------------------------------------
221      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
222      !!
223      INTEGER ::   jn   ! dummy loop index
224      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
225      !!
226      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,   &
227         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,   &
228         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx  , sn_rcv_co2
229      !!---------------------------------------------------------------------
230      !
231      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
232      !
233      CALL wrk_alloc( jpi,jpj, zacs, zaos )
234
235      ! ================================ !
236      !      Namelist informations       !
237      ! ================================ !
238
239      ! default definitions
240      !                    !     description       !  multiple  !    vector   !      vector          ! vector !
241      !                    !                       ! categories !  reference  !    orientation       ! grids  !
242      ! send
243      sn_snd_temp   = FLD_C( 'weighted oce and ice',    'no'    ,     ''      ,         ''           ,   ''   ) 
244      sn_snd_alb    = FLD_C( 'weighted ice'        ,    'no'    ,     ''      ,         ''           ,   ''   ) 
245      sn_snd_thick  = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''           ,   ''   ) 
246      sn_snd_crt    = FLD_C( 'none'                ,    'no'    , 'spherical' , 'eastward-northward' ,  'T'   )     
247      sn_snd_co2    = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''           ,   ''   )     
248      ! receive
249      sn_rcv_w10m   = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
250      sn_rcv_taumod = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
251      sn_rcv_tau    = FLD_C( 'oce only'            ,    'no'    , 'cartesian' , 'eastward-northward',  'U,V'  ) 
252      sn_rcv_dqnsdt = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
253      sn_rcv_qsr    = FLD_C( 'oce and ice'         ,    'no'    ,     ''      ,         ''          ,   ''    )
254      sn_rcv_qns    = FLD_C( 'oce and ice'         ,    'no'    ,     ''      ,         ''          ,   ''    )
255      sn_rcv_emp    = FLD_C( 'conservative'        ,    'no'    ,     ''      ,         ''          ,   ''    )
256      sn_rcv_rnf    = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
257      sn_rcv_cal    = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
258      sn_rcv_iceflx = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
259      sn_rcv_co2    = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
260
261      REWIND( numnam )                    ! ... read namlist namsbc_cpl
262      READ  ( numnam, namsbc_cpl )
263
264      IF(lwp) THEN                        ! control print
265         WRITE(numout,*)
266         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
267         WRITE(numout,*)'~~~~~~~~~~~~'
268         WRITE(numout,*)'  received fields (mutiple ice categogies)'
269         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
270         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
271         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
272         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
273         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
274         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
275         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
276         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
277         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
278         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
279         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
280         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
281         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
282         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
283         WRITE(numout,*)'  sent fields (multiple ice categories)'
284         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
285         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
286         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
287         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
288         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
289         WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
290         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
291         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
292      ENDIF
293
294      !                                   ! allocate sbccpl arrays
295      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
296     
297      ! ================================ !
298      !   Define the receive interface   !
299      ! ================================ !
300      nrcvinfo(:) = OASIS_idle   ! needed by nrcvinfo(jpr_otx1) if we do not receive ocean stress
301
302      ! for each field: define the OASIS name                              (srcv(:)%clname)
303      !                 define receive or not from the namelist parameters (srcv(:)%laction)
304      !                 define the north fold type of lbc                  (srcv(:)%nsgn)
305
306      ! default definitions of srcv
307      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
308
309      !                                                      ! ------------------------- !
310      !                                                      ! ice and ocean wind stress !   
311      !                                                      ! ------------------------- !
312      !                                                           ! Name
313      srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
314      srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
315      srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
316      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
317      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
318      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
319      !
320      srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
321      srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
322      srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
323      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
324      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
325      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
326      !
327      ! Vectors: change of sign at north fold ONLY if on the local grid
328      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
329     
330      !                                                           ! Set grid and action
331      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
332      CASE( 'T' ) 
333         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
334         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
335         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
336      CASE( 'U,V' ) 
337         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
338         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
339         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
340         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
341         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
342      CASE( 'U,V,T' )
343         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
344         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
345         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
346         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
347         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
348      CASE( 'U,V,I' )
349         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
350         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
351         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
352         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
353         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
354      CASE( 'U,V,F' )
355         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
356         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
357         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
358         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
359         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
360      CASE( 'T,I' ) 
361         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
362         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
363         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
364         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
365      CASE( 'T,F' ) 
366         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
367         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
368         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
369         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
370      CASE( 'T,U,V' )
371         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
372         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
373         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
374         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
375         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
376      CASE default   
377         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
378      END SELECT
379      !
380      IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
381         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE. 
382      !
383      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
384         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
385         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
386         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
387      ENDIF
388       
389      !                                                      ! ------------------------- !
390      !                                                      !    freshwater budget      !   E-P
391      !                                                      ! ------------------------- !
392      ! we suppose that atmosphere modele do not make the difference between precipiration (liquide or solid)
393      ! over ice of free ocean within the same atmospheric cell.cd
394      srcv(jpr_rain)%clname = 'OTotRain'      ! Rain = liquid precipitation
395      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
396      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
397      srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation
398      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation
399      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
400      srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
401      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
402      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
403      CASE( 'conservative'  )   ;   srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
404      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
405      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
406      END SELECT
407
408      !                                                      ! ------------------------- !
409      !                                                      !     Runoffs & Calving     !   
410      !                                                      ! ------------------------- !
411      srcv(jpr_rnf   )%clname = 'O_Runoff'   ;   IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' )   srcv(jpr_rnf)%laction = .TRUE.
412! This isn't right - really just want ln_rnf_emp changed
413!                                                 IF( TRIM( sn_rcv_rnf%cldes ) == 'climato' )   THEN   ;   ln_rnf = .TRUE.
414!                                                 ELSE                                                 ;   ln_rnf = .FALSE.
415!                                                 ENDIF
416      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
417
418      !                                                      ! ------------------------- !
419      !                                                      !    non solar radiation    !   Qns
420      !                                                      ! ------------------------- !
421      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
422      srcv(jpr_qnsice)%clname = 'O_QnsIce'
423      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
424      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
425      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
426      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
427      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
428      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
429      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
430      END SELECT
431      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
432         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
433      !                                                      ! ------------------------- !
434      !                                                      !    solar radiation        !   Qsr
435      !                                                      ! ------------------------- !
436      srcv(jpr_qsroce)%clname = 'O_QsrOce'
437      srcv(jpr_qsrice)%clname = 'O_QsrIce'
438      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
439      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
440      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
441      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
442      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
443      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
444      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
445      END SELECT
446      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
447         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
448      !                                                      ! ------------------------- !
449      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
450      !                                                      ! ------------------------- !
451      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
452      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
453      !
454      ! non solar sensitivity mandatory for LIM ice model
455      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4) &
456         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
457      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
458      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
459         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
460      !                                                      ! ------------------------- !
461      !                                                      !    Ice Qsr penetration    !   
462      !                                                      ! ------------------------- !
463      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer
464      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
465      ! Coupled case: since cloud cover is not received from atmosphere
466      !               ===> defined as constant value -> definition done in sbc_cpl_init
467      fr1_i0(:,:) = 0.18
468      fr2_i0(:,:) = 0.82
469      !                                                      ! ------------------------- !
470      !                                                      !      10m wind module      !   
471      !                                                      ! ------------------------- !
472      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
473      !
474      !                                                      ! ------------------------- !
475      !                                                      !   wind stress module      !   
476      !                                                      ! ------------------------- !
477      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
478      lhftau = srcv(jpr_taum)%laction
479
480      !                                                      ! ------------------------- !
481      !                                                      !      Atmospheric CO2      !
482      !                                                      ! ------------------------- !
483      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
484      !                                                      ! ------------------------- !
485      !                                                      !   topmelt and botmelt     !   
486      !                                                      ! ------------------------- !
487      srcv(jpr_topm )%clname = 'OTopMlt'
488      srcv(jpr_botm )%clname = 'OBotMlt'
489      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
490         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
491            srcv(jpr_topm:jpr_botm)%nct = jpl
492         ELSE
493            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
494         ENDIF
495         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
496      ENDIF
497
498      ! Allocate all parts of frcv used for received fields
499      DO jn = 1, jprcv
500         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
501      END DO
502      ! Allocate taum part of frcv which is used even when not received as coupling field
503      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jn)%nct) )
504
505      ! ================================ !
506      !     Define the send interface    !
507      ! ================================ !
508      ! for each field: define the OASIS name                           (ssnd(:)%clname)
509      !                 define send or not from the namelist parameters (ssnd(:)%laction)
510      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
511     
512      ! default definitions of nsnd
513      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
514         
515      !                                                      ! ------------------------- !
516      !                                                      !    Surface temperature    !
517      !                                                      ! ------------------------- !
518      ssnd(jps_toce)%clname = 'O_SSTSST'
519      ssnd(jps_tice)%clname = 'O_TepIce'
520      ssnd(jps_tmix)%clname = 'O_TepMix'
521      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
522      CASE( 'oce only'             )   ;   ssnd(   jps_toce             )%laction = .TRUE.
523      CASE( 'weighted oce and ice' )
524         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
525         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
526      CASE( 'mixed oce-ice'        )   ;   ssnd(   jps_tmix             )%laction = .TRUE.
527      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
528      END SELECT
529     
530      !                                                      ! ------------------------- !
531      !                                                      !          Albedo           !
532      !                                                      ! ------------------------- !
533      ssnd(jps_albice)%clname = 'O_AlbIce' 
534      ssnd(jps_albmix)%clname = 'O_AlbMix'
535      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
536      CASE( 'none'          )       ! nothing to do
537      CASE( 'weighted ice'  )   ;   ssnd(jps_albice)%laction = .TRUE.
538      CASE( 'mixed oce-ice' )   ;   ssnd(jps_albmix)%laction = .TRUE.
539      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
540      END SELECT
541      !
542      ! Need to calculate oceanic albedo if
543      !     1. sending mixed oce-ice albedo or
544      !     2. receiving mixed oce-ice solar radiation
545      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
546         CALL albedo_oce( zaos, zacs )
547         ! Due to lack of information on nebulosity : mean clear/overcast sky
548         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
549      ENDIF
550
551      !                                                      ! ------------------------- !
552      !                                                      !  Ice fraction & Thickness !
553      !                                                      ! ------------------------- !
554      ssnd(jps_fice)%clname = 'OIceFrc'
555      ssnd(jps_hice)%clname = 'OIceTck'
556      ssnd(jps_hsnw)%clname = 'OSnwTck'
557      IF( k_ice /= 0 ) THEN
558         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
559! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
560         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
561      ENDIF
562
563      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
564      CASE ( 'ice and snow' ) 
565         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
566         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
567            ssnd(jps_hice:jps_hsnw)%nct = jpl
568         ELSE
569            IF ( jpl > 1 ) THEN
570               CALL ctl_stop( 'sbc_cpl_init: use weighted ice and snow option for sn_snd_thick%cldes if not exchanging category fields' )
571            ENDIF
572         ENDIF
573      CASE ( 'weighted ice and snow' ) 
574         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
575         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
576      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
577      END SELECT
578
579      !                                                      ! ------------------------- !
580      !                                                      !      Surface current      !
581      !                                                      ! ------------------------- !
582      !        ocean currents              !            ice velocities
583      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
584      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
585      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
586      !
587      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
588
589      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
590         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
591      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
592         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
593         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
594      ENDIF
595      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
596      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
597      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
598      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
599      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
600      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
601      CASE( 'weighted oce and ice' )   !   nothing to do
602      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
603      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
604      END SELECT
605
606      !                                                      ! ------------------------- !
607      !                                                      !          CO2 flux         !
608      !                                                      ! ------------------------- !
609      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
610      !
611      ! ================================ !
612      !   initialisation of the coupler  !
613      ! ================================ !
614
615      CALL cpl_prism_define(jprcv, jpsnd)           
616      !
617      IF( ln_dm2dc .AND. ( cpl_prism_freq( jpr_qsroce ) + cpl_prism_freq( jpr_qsrmix ) /= 86400 ) )   &
618         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
619
620      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
621      !
622      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
623      !
624   END SUBROUTINE sbc_cpl_init
625
626
627   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
628      !!----------------------------------------------------------------------
629      !!             ***  ROUTINE sbc_cpl_rcv  ***
630      !!
631      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
632      !!                provide the ocean heat and freshwater fluxes.
633      !!
634      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
635      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
636      !!                to know if the field was really received or not
637      !!
638      !!              --> If ocean stress was really received:
639      !!
640      !!                  - transform the received ocean stress vector from the received
641      !!                 referential and grid into an atmosphere-ocean stress in
642      !!                 the (i,j) ocean referencial and at the ocean velocity point.
643      !!                    The received stress are :
644      !!                     - defined by 3 components (if cartesian coordinate)
645      !!                            or by 2 components (if spherical)
646      !!                     - oriented along geographical   coordinate (if eastward-northward)
647      !!                            or  along the local grid coordinate (if local grid)
648      !!                     - given at U- and V-point, resp.   if received on 2 grids
649      !!                            or at T-point               if received on 1 grid
650      !!                    Therefore and if necessary, they are successively
651      !!                  processed in order to obtain them
652      !!                     first  as  2 components on the sphere
653      !!                     second as  2 components oriented along the local grid
654      !!                     third  as  2 components on the U,V grid
655      !!
656      !!              -->
657      !!
658      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
659      !!             and total ocean freshwater fluxes 
660      !!
661      !! ** Method  :   receive all fields from the atmosphere and transform
662      !!              them into ocean surface boundary condition fields
663      !!
664      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
665      !!                        taum, wndm   wind stres and wind speed module at T-point
666      !!                        qns , qsr    non solar and solar ocean heat fluxes   ('ocean only case)
667      !!                        emp = emps   evap. - precip. (- runoffs) (- calving) ('ocean only case)
668      !!----------------------------------------------------------------------
669      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
670      INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation
671      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
672      !!
673      LOGICAL ::    llnewtx, llnewtau      ! update wind stress components and module??
674      INTEGER  ::   ji, jj, jn             ! dummy loop indices
675      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
676      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
677      REAL(wp) ::   zcoef                  ! temporary scalar
678      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
679      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
680      REAL(wp) ::   zzx, zzy               ! temporary variables
681      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
682      !!----------------------------------------------------------------------
683      !
684      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
685      !
686      CALL wrk_alloc( jpi,jpj, ztx, zty )
687
688      IF( kt == nit000 )   CALL sbc_cpl_init( k_ice )          ! initialisation
689
690      !                                                 ! Receive all the atmos. fields (including ice information)
691      isec = ( kt - nit000 ) * NINT( rdttra(1) )             ! date of exchanges
692      DO jn = 1, jprcv                                       ! received fields sent by the atmosphere
693         IF( srcv(jn)%laction )   CALL cpl_prism_rcv( jn, isec, frcv(jn)%z3, nrcvinfo(jn) )
694      END DO
695
696      !                                                      ! ========================= !
697      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
698         !                                                   ! ========================= !
699         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
700         ! => need to be done only when we receive the field
701         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
702            !
703            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
704               !                                                       ! (cartesian to spherical -> 3 to 2 components)
705               !
706               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
707                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
708               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
709               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
710               !
711               IF( srcv(jpr_otx2)%laction ) THEN
712                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
713                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
714                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
715                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
716               ENDIF
717               !
718            ENDIF
719            !
720            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
721               !                                                       ! (geographical to local grid -> rotate the components)
722               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
723               IF( srcv(jpr_otx2)%laction ) THEN
724                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
725               ELSE 
726                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
727               ENDIF
728               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
729               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
730            ENDIF
731            !                             
732            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
733               DO jj = 2, jpjm1                                          ! T ==> (U,V)
734                  DO ji = fs_2, fs_jpim1   ! vector opt.
735                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
736                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
737                  END DO
738               END DO
739               CALL lbc_lnk( frcv(jpr_otx1)%z3(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
740            ENDIF
741            llnewtx = .TRUE.
742         ELSE
743            llnewtx = .FALSE.
744         ENDIF
745         !                                                   ! ========================= !
746      ELSE                                                   !   No dynamical coupling   !
747         !                                                   ! ========================= !
748         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
749         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
750         llnewtx = .TRUE.
751         !
752      ENDIF
753     
754      !                                                      ! ========================= !
755      !                                                      !    wind stress module     !   (taum)
756      !                                                      ! ========================= !
757      !
758      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
759         ! => need to be done only when otx1 was changed
760         IF( llnewtx ) THEN
761!CDIR NOVERRCHK
762            DO jj = 2, jpjm1
763!CDIR NOVERRCHK
764               DO ji = fs_2, fs_jpim1   ! vect. opt.
765                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
766                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
767                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
768               END DO
769            END DO
770            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
771            llnewtau = .TRUE.
772         ELSE
773            llnewtau = .FALSE.
774         ENDIF
775      ELSE
776         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
777         ! Stress module can be negative when received (interpolation problem)
778         IF( llnewtau ) THEN
779            frcv(jpr_taum)%z3(:,:,1) = MAX( 0.0e0, frcv(jpr_taum)%z3(:,:,1) )
780         ENDIF
781      ENDIF
782     
783      !                                                      ! ========================= !
784      !                                                      !      10 m wind speed      !   (wndm)
785      !                                                      ! ========================= !
786      !
787      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
788         ! => need to be done only when taumod was changed
789         IF( llnewtau ) THEN
790            zcoef = 1. / ( zrhoa * zcdrag ) 
791!CDIR NOVERRCHK
792            DO jj = 1, jpj
793!CDIR NOVERRCHK
794               DO ji = 1, jpi 
795                  wndm(ji,jj) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
796               END DO
797            END DO
798         ENDIF
799      ELSE
800         IF ( nrcvinfo(jpr_w10m) == OASIS_Rcv ) wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
801      ENDIF
802
803      ! u(v)tau and taum will be modified by ice model
804      ! -> need to be reset before each call of the ice/fsbc     
805      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
806         !
807         utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
808         vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
809         taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
810         CALL iom_put( "taum_oce", taum )   ! output wind stress module
811         
812      ENDIF
813
814#if defined key_cpl_carbon_cycle
815      !                                                              ! atmosph. CO2 (ppm)
816      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
817#endif
818
819      !                                                      ! ========================= !
820      IF( k_ice <= 1 ) THEN                                  !  heat & freshwater fluxes ! (Ocean only case)
821         !                                                   ! ========================= !
822         !
823         !                                                       ! non solar heat flux over the ocean (qns)
824         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
825         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
826         ! add the latent heat of solid precip. melting
827         IF( srcv(jpr_snow  )%laction )   qns(:,:) = qns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus             
828
829         !                                                       ! solar flux over the ocean          (qsr)
830         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
831         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
832         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
833         !
834         !                                                       ! total freshwater fluxes over the ocean (emp, emps)
835         SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
836         CASE( 'conservative' )
837            emp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
838         CASE( 'oce only', 'oce and ice' )
839            emp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
840         CASE default
841            CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
842         END SELECT
843         !
844         !                                                        ! runoffs and calving (added in emp)
845         IF( srcv(jpr_rnf)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
846         IF( srcv(jpr_cal)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_cal)%z3(:,:,1)
847         !
848!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
849!!gm                                       at least should be optional...
850!!         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN     ! add to the total freshwater budget
851!!            ! remove negative runoff
852!!            zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
853!!            zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
854!!            IF( lk_mpp )   CALL mpp_sum( zcumulpos )   ! sum over the global domain
855!!            IF( lk_mpp )   CALL mpp_sum( zcumulneg )
856!!            IF( zcumulpos /= 0. ) THEN                 ! distribute negative runoff on positive runoff grid points
857!!               zcumulneg = 1.e0 + zcumulneg / zcumulpos
858!!               frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
859!!            ENDIF     
860!!            ! add runoff to e-p
861!!            emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
862!!         ENDIF
863!!gm  end of internal cooking
864         !
865         emps(:,:) = emp(:,:)                                        ! concentration/dilution = emp
866 
867      ENDIF
868      !
869      CALL wrk_dealloc( jpi,jpj, ztx, zty )
870      !
871      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
872      !
873   END SUBROUTINE sbc_cpl_rcv
874   
875
876   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
877      !!----------------------------------------------------------------------
878      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
879      !!
880      !! ** Purpose :   provide the stress over sea-ice in coupled mode
881      !!
882      !! ** Method  :   transform the received stress from the atmosphere into
883      !!             an atmosphere-ice stress in the (i,j) ocean referencial
884      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
885      !!                'C'-grid : i- (j-) components given at U- (V-) point
886      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
887      !!
888      !!                The received stress are :
889      !!                 - defined by 3 components (if cartesian coordinate)
890      !!                        or by 2 components (if spherical)
891      !!                 - oriented along geographical   coordinate (if eastward-northward)
892      !!                        or  along the local grid coordinate (if local grid)
893      !!                 - given at U- and V-point, resp.   if received on 2 grids
894      !!                        or at a same point (T or I) if received on 1 grid
895      !!                Therefore and if necessary, they are successively
896      !!             processed in order to obtain them
897      !!                 first  as  2 components on the sphere
898      !!                 second as  2 components oriented along the local grid
899      !!                 third  as  2 components on the cp_ice_msh point
900      !!
901      !!                In 'oce and ice' case, only one vector stress field
902      !!             is received. It has already been processed in sbc_cpl_rcv
903      !!             so that it is now defined as (i,j) components given at U-
904      !!             and V-points, respectively. Therefore, here only the third
905      !!             transformation is done and only if the ice-grid is a 'I'-grid.
906      !!
907      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
908      !!----------------------------------------------------------------------
909      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
910      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
911      !!
912      INTEGER ::   ji, jj                          ! dummy loop indices
913      INTEGER ::   itx                             ! index of taux over ice
914      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
915      !!----------------------------------------------------------------------
916      !
917      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
918      !
919      CALL wrk_alloc( jpi,jpj, ztx, zty )
920
921      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
922      ELSE                                ;   itx =  jpr_otx1
923      ENDIF
924
925      ! do something only if we just received the stress from atmosphere
926      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
927
928         !                                                      ! ======================= !
929         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
930            !                                                   ! ======================= !
931           
932            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
933               !                                                       ! (cartesian to spherical -> 3 to 2 components)
934               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
935                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
936               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
937               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
938               !
939               IF( srcv(jpr_itx2)%laction ) THEN
940                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
941                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
942                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
943                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
944               ENDIF
945               !
946            ENDIF
947            !
948            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
949               !                                                       ! (geographical to local grid -> rotate the components)
950               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
951               IF( srcv(jpr_itx2)%laction ) THEN
952                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
953               ELSE
954                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
955               ENDIF
956               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
957               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
958            ENDIF
959            !                                                   ! ======================= !
960         ELSE                                                   !     use ocean stress    !
961            !                                                   ! ======================= !
962            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
963            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
964            !
965         ENDIF
966
967         !                                                      ! ======================= !
968         !                                                      !     put on ice grid     !
969         !                                                      ! ======================= !
970         !   
971         !                                                  j+1   j     -----V---F
972         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
973         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
974         !                                                               |       |
975         !                                                   j    j-1   -I-------|
976         !                                               (for I)         |       |
977         !                                                              i-1  i   i
978         !                                                               i      i+1 (for I)
979         SELECT CASE ( cp_ice_msh )
980            !
981         CASE( 'I' )                                         ! B-grid ==> I
982            SELECT CASE ( srcv(jpr_itx1)%clgrid )
983            CASE( 'U' )
984               DO jj = 2, jpjm1                                   ! (U,V) ==> I
985                  DO ji = 2, jpim1   ! NO vector opt.
986                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
987                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
988                  END DO
989               END DO
990            CASE( 'F' )
991               DO jj = 2, jpjm1                                   ! F ==> I
992                  DO ji = 2, jpim1   ! NO vector opt.
993                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
994                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
995                  END DO
996               END DO
997            CASE( 'T' )
998               DO jj = 2, jpjm1                                   ! T ==> I
999                  DO ji = 2, jpim1   ! NO vector opt.
1000                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
1001                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
1002                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
1003                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1004                  END DO
1005               END DO
1006            CASE( 'I' )
1007               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
1008               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1009            END SELECT
1010            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1011               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1012            ENDIF
1013            !
1014         CASE( 'F' )                                         ! B-grid ==> F
1015            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1016            CASE( 'U' )
1017               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1018                  DO ji = fs_2, fs_jpim1   ! vector opt.
1019                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1020                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1021                  END DO
1022               END DO
1023            CASE( 'I' )
1024               DO jj = 2, jpjm1                                   ! I ==> F
1025                  DO ji = 2, jpim1   ! NO vector opt.
1026                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1027                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1028                  END DO
1029               END DO
1030            CASE( 'T' )
1031               DO jj = 2, jpjm1                                   ! T ==> F
1032                  DO ji = 2, jpim1   ! NO vector opt.
1033                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1034                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1035                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1036                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1037                  END DO
1038               END DO
1039            CASE( 'F' )
1040               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1041               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1042            END SELECT
1043            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1044               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1045            ENDIF
1046            !
1047         CASE( 'C' )                                         ! C-grid ==> U,V
1048            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1049            CASE( 'U' )
1050               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1051               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1052            CASE( 'F' )
1053               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1054                  DO ji = fs_2, fs_jpim1   ! vector opt.
1055                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1056                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1057                  END DO
1058               END DO
1059            CASE( 'T' )
1060               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1061                  DO ji = fs_2, fs_jpim1   ! vector opt.
1062                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1063                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1064                  END DO
1065               END DO
1066            CASE( 'I' )
1067               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1068                  DO ji = 2, jpim1   ! NO vector opt.
1069                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1070                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1071                  END DO
1072               END DO
1073            END SELECT
1074            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1075               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1076            ENDIF
1077         END SELECT
1078
1079      ENDIF
1080      !   
1081      CALL wrk_dealloc( jpi,jpj, ztx, zty )
1082      !
1083      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
1084      !
1085   END SUBROUTINE sbc_cpl_ice_tau
1086   
1087
1088   SUBROUTINE sbc_cpl_ice_flx( p_frld  , palbi   , psst    , pist    )
1089      !!----------------------------------------------------------------------
1090      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1091      !!
1092      !! ** Purpose :   provide the heat and freshwater fluxes of the
1093      !!              ocean-ice system.
1094      !!
1095      !! ** Method  :   transform the fields received from the atmosphere into
1096      !!             surface heat and fresh water boundary condition for the
1097      !!             ice-ocean system. The following fields are provided:
1098      !!              * total non solar, solar and freshwater fluxes (qns_tot,
1099      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1100      !!             NB: emp_tot include runoffs and calving.
1101      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1102      !!             emp_ice = sublimation - solid precipitation as liquid
1103      !!             precipitation are re-routed directly to the ocean and
1104      !!             runoffs and calving directly enter the ocean.
1105      !!              * solid precipitation (sprecip), used to add to qns_tot
1106      !!             the heat lost associated to melting solid precipitation
1107      !!             over the ocean fraction.
1108      !!       ===>> CAUTION here this changes the net heat flux received from
1109      !!             the atmosphere
1110      !!
1111      !!                  - the fluxes have been separated from the stress as
1112      !!                 (a) they are updated at each ice time step compare to
1113      !!                 an update at each coupled time step for the stress, and
1114      !!                 (b) the conservative computation of the fluxes over the
1115      !!                 sea-ice area requires the knowledge of the ice fraction
1116      !!                 after the ice advection and before the ice thermodynamics,
1117      !!                 so that the stress is updated before the ice dynamics
1118      !!                 while the fluxes are updated after it.
1119      !!
1120      !! ** Action  :   update at each nf_ice time step:
1121      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1122      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1123      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
1124      !!                   emp_ice            ice sublimation - solid precipitation over the ice
1125      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
1126      !!                   sprecip             solid precipitation over the ocean 
1127      !!----------------------------------------------------------------------
1128      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1129      ! optional arguments, used only in 'mixed oce-ice' case
1130      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1131      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celcius]
1132      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature     [Kelvin]
1133      !
1134      INTEGER ::   jl   ! dummy loop index
1135      REAL(wp), POINTER, DIMENSION(:,:) ::   zcptn, ztmp, zicefr
1136      !!----------------------------------------------------------------------
1137      !
1138      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1139      !
1140      CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr )
1141
1142      zicefr(:,:) = 1.- p_frld(:,:)
1143      IF( lk_diaar5 )   zcptn(:,:) = rcp * tsn(:,:,1,jp_tem)
1144      !
1145      !                                                      ! ========================= !
1146      !                                                      !    freshwater budget      !   (emp)
1147      !                                                      ! ========================= !
1148      !
1149      !                                                           ! total Precipitations - total Evaporation (emp_tot)
1150      !                                                           ! solid precipitation  - sublimation       (emp_ice)
1151      !                                                           ! solid Precipitation                      (sprecip)
1152      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1153      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1154         sprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                 ! May need to ensure positive here
1155         tprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + sprecip (:,:) ! May need to ensure positive here
1156         emp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - tprecip(:,:)
1157         emp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1158                           CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation
1159         IF( lk_diaar5 )   CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
1160         ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1161                           CALL iom_put( 'evap_ao_cea'  , ztmp                            )   ! ice-free oce evap (cell average)
1162         IF( lk_diaar5 )   CALL iom_put( 'hflx_evap_cea', ztmp(:,:         ) * zcptn(:,:) )   ! heat flux from from evap (cell ave)
1163      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1164         emp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1165         emp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
1166         sprecip(:,:) = - frcv(jpr_semp)%z3(:,:,1) + frcv(jpr_ievp)%z3(:,:,1)
1167      END SELECT
1168
1169      CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
1170      CALL iom_put( 'snow_ao_cea', sprecip(:,:         ) * p_frld(:,:)    )   ! Snow        over ice-free ocean  (cell average)
1171      CALL iom_put( 'snow_ai_cea', sprecip(:,:         ) * zicefr(:,:)    )   ! Snow        over sea-ice         (cell average)
1172      CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
1173      !   
1174      !                                                           ! runoffs and calving (put in emp_tot)
1175      IF( srcv(jpr_rnf)%laction ) THEN
1176         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)
1177                           CALL iom_put( 'runoffs'      , frcv(jpr_rnf)%z3(:,:,1)              )   ! rivers
1178         IF( lk_diaar5 )   CALL iom_put( 'hflx_rnf_cea' , frcv(jpr_rnf)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from rivers
1179      ENDIF
1180      IF( srcv(jpr_cal)%laction ) THEN
1181         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1182         CALL iom_put( 'calving', frcv(jpr_cal)%z3(:,:,1) )
1183      ENDIF
1184      !
1185!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
1186!!gm                                       at least should be optional...
1187!!       ! remove negative runoff                            ! sum over the global domain
1188!!       zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
1189!!       zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
1190!!       IF( lk_mpp )   CALL mpp_sum( zcumulpos )
1191!!       IF( lk_mpp )   CALL mpp_sum( zcumulneg )
1192!!       IF( zcumulpos /= 0. ) THEN                          ! distribute negative runoff on positive runoff grid points
1193!!          zcumulneg = 1.e0 + zcumulneg / zcumulpos
1194!!          frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
1195!!       ENDIF     
1196!!       emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)   ! add runoff to e-p
1197!!
1198!!gm  end of internal cooking
1199
1200      !                                                      ! ========================= !
1201      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1202      !                                                      ! ========================= !
1203      CASE( 'oce only' )                                     ! the required field is directly provided
1204         qns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
1205      CASE( 'conservative' )                                      ! the required fields are directly provided
1206         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1207         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1208            qns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1209         ELSE
1210            ! Set all category values equal for the moment
1211            DO jl=1,jpl
1212               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1213            ENDDO
1214         ENDIF
1215      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
1216         qns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1217         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1218            DO jl=1,jpl
1219               qns_tot(:,:   ) = qns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1220               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1221            ENDDO
1222         ELSE
1223            DO jl=1,jpl
1224               qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1225               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1226            ENDDO
1227         ENDIF
1228      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
1229! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1230         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1231         qns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1232            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1233            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
1234      END SELECT
1235      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus               ! add the latent heat of solid precip. melting
1236      qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)                     ! over free ocean
1237      IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1238!!gm
1239!!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in
1240!!    the flux that enter the ocean....
1241!!    moreover 1 - it is not diagnose anywhere....
1242!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1243!!
1244!! similar job should be done for snow and precipitation temperature
1245      !                                     
1246      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting
1247         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting
1248         qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)
1249         IF( lk_diaar5 )   CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1250      ENDIF
1251
1252      !                                                      ! ========================= !
1253      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1254      !                                                      ! ========================= !
1255      CASE( 'oce only' )
1256         qsr_tot(:,:  ) = MAX(0.0,frcv(jpr_qsroce)%z3(:,:,1))
1257      CASE( 'conservative' )
1258         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1259         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1260            qsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1261         ELSE
1262            ! Set all category values equal for the moment
1263            DO jl=1,jpl
1264               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1265            ENDDO
1266         ENDIF
1267         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1268         qsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1269      CASE( 'oce and ice' )
1270         qsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1271         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1272            DO jl=1,jpl
1273               qsr_tot(:,:   ) = qsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1274               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1275            ENDDO
1276         ELSE
1277            DO jl=1,jpl
1278               qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1279               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1280            ENDDO
1281         ENDIF
1282      CASE( 'mixed oce-ice' )
1283         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1284! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1285!       Create solar heat flux over ice using incoming solar heat flux and albedos
1286!       ( see OASIS3 user guide, 5th edition, p39 )
1287         qsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1288            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1289            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1290      END SELECT
1291      IF( ln_dm2dc ) THEN   ! modify qsr to include the diurnal cycle
1292         qsr_tot(:,:  ) = sbc_dcy( qsr_tot(:,:  ) )
1293         DO jl=1,jpl
1294            qsr_ice(:,:,jl) = sbc_dcy( qsr_ice(:,:,jl) )
1295         ENDDO
1296      ENDIF
1297
1298      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )
1299      CASE ('coupled')
1300         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1301            dqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
1302         ELSE
1303            ! Set all category values equal for the moment
1304            DO jl=1,jpl
1305               dqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
1306            ENDDO
1307         ENDIF
1308      END SELECT
1309
1310      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )
1311      CASE ('coupled')
1312         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
1313         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
1314      END SELECT
1315
1316      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
1317      !
1318      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
1319      !
1320   END SUBROUTINE sbc_cpl_ice_flx
1321   
1322   
1323   SUBROUTINE sbc_cpl_snd( kt )
1324      !!----------------------------------------------------------------------
1325      !!             ***  ROUTINE sbc_cpl_snd  ***
1326      !!
1327      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
1328      !!
1329      !! ** Method  :   send to the atmosphere through a call to cpl_prism_snd
1330      !!              all the needed fields (as defined in sbc_cpl_init)
1331      !!----------------------------------------------------------------------
1332      INTEGER, INTENT(in) ::   kt
1333      !
1334      INTEGER ::   ji, jj, jl   ! dummy loop indices
1335      INTEGER ::   isec, info   ! local integer
1336      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
1337      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
1338      !!----------------------------------------------------------------------
1339      !
1340      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
1341      !
1342      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1343      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1344
1345      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
1346
1347      zfr_l(:,:) = 1.- fr_i(:,:)
1348
1349      !                                                      ! ------------------------- !
1350      !                                                      !    Surface temperature    !   in Kelvin
1351      !                                                      ! ------------------------- !
1352      SELECT CASE( sn_snd_temp%cldes)
1353      CASE( 'oce only'             )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0
1354      CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:)   
1355         SELECT CASE( sn_snd_temp%clcat )
1356         CASE( 'yes' )   
1357            ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1358         CASE( 'no' )
1359            ztmp3(:,:,:) = 0.0
1360            DO jl=1,jpl
1361               ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
1362            ENDDO
1363         CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1364         END SELECT
1365      CASE( 'mixed oce-ice'        )   
1366         ztmp1(:,:) = ( tsn(:,:,1,1) + rt0 ) * zfr_l(:,:) 
1367         DO jl=1,jpl
1368            ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
1369         ENDDO
1370      CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
1371      END SELECT
1372      IF( ssnd(jps_toce)%laction )   CALL cpl_prism_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1373      IF( ssnd(jps_tice)%laction )   CALL cpl_prism_snd( jps_tice, isec, ztmp3, info )
1374      IF( ssnd(jps_tmix)%laction )   CALL cpl_prism_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1375      !
1376      !                                                      ! ------------------------- !
1377      !                                                      !           Albedo          !
1378      !                                                      ! ------------------------- !
1379      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
1380         ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1381         CALL cpl_prism_snd( jps_albice, isec, ztmp3, info )
1382      ENDIF
1383      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
1384         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
1385         DO jl=1,jpl
1386            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
1387         ENDDO
1388         CALL cpl_prism_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1389      ENDIF
1390      !                                                      ! ------------------------- !
1391      !                                                      !  Ice fraction & Thickness !
1392      !                                                      ! ------------------------- !
1393      ! Send ice fraction field
1394      SELECT CASE( sn_snd_thick%clcat )
1395         CASE( 'yes' )   
1396            ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
1397         CASE( 'no' )
1398            ztmp3(:,:,1) = fr_i(:,:)
1399      CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1400      END SELECT
1401      IF( ssnd(jps_fice)%laction ) CALL cpl_prism_snd( jps_fice, isec, ztmp3, info )
1402
1403      ! Send ice and snow thickness field
1404      SELECT CASE( sn_snd_thick%cldes)
1405      CASE( 'weighted ice and snow' )   
1406         SELECT CASE( sn_snd_thick%clcat )
1407         CASE( 'yes' )   
1408            ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
1409            ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
1410         CASE( 'no' )
1411            ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
1412            DO jl=1,jpl
1413               ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
1414               ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
1415            ENDDO
1416         CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1417         END SELECT
1418      CASE( 'ice and snow'         )   
1419         ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
1420         ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
1421      CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
1422      END SELECT
1423      IF( ssnd(jps_hice)%laction )   CALL cpl_prism_snd( jps_hice, isec, ztmp3, info )
1424      IF( ssnd(jps_hsnw)%laction )   CALL cpl_prism_snd( jps_hsnw, isec, ztmp4, info )
1425      !
1426#if defined key_cpl_carbon_cycle
1427      !                                                      ! ------------------------- !
1428      !                                                      !  CO2 flux from PISCES     !
1429      !                                                      ! ------------------------- !
1430      IF( ssnd(jps_co2)%laction )   CALL cpl_prism_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
1431      !
1432#endif
1433      !                                                      ! ------------------------- !
1434      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
1435         !                                                   ! ------------------------- !
1436         !   
1437         !                                                  j+1   j     -----V---F
1438         ! surface velocity always sent from T point                     !       |
1439         !                                                        j      |   T   U
1440         !                                                               |       |
1441         !                                                   j    j-1   -I-------|
1442         !                                               (for I)         |       |
1443         !                                                              i-1  i   i
1444         !                                                               i      i+1 (for I)
1445         SELECT CASE( TRIM( sn_snd_crt%cldes ) )
1446         CASE( 'oce only'             )      ! C-grid ==> T
1447            DO jj = 2, jpjm1
1448               DO ji = fs_2, fs_jpim1   ! vector opt.
1449                  zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
1450                  zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
1451               END DO
1452            END DO
1453         CASE( 'weighted oce and ice' )   
1454            SELECT CASE ( cp_ice_msh )
1455            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1456               DO jj = 2, jpjm1
1457                  DO ji = fs_2, fs_jpim1   ! vector opt.
1458                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1459                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
1460                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1461                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1462                  END DO
1463               END DO
1464            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1465               DO jj = 2, jpjm1
1466                  DO ji = 2, jpim1   ! NO vector opt.
1467                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1468                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1469                     zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1470                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1471                     zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1472                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1473                  END DO
1474               END DO
1475            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1476               DO jj = 2, jpjm1
1477                  DO ji = 2, jpim1   ! NO vector opt.
1478                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1479                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1480                     zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1481                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1482                     zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1483                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1484                  END DO
1485               END DO
1486            END SELECT
1487            CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
1488         CASE( 'mixed oce-ice'        )
1489            SELECT CASE ( cp_ice_msh )
1490            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1491               DO jj = 2, jpjm1
1492                  DO ji = fs_2, fs_jpim1   ! vector opt.
1493                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
1494                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1495                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
1496                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1497                  END DO
1498               END DO
1499            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1500               DO jj = 2, jpjm1
1501                  DO ji = 2, jpim1   ! NO vector opt.
1502                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1503                        &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1504                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1505                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1506                        &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1507                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1508                  END DO
1509               END DO
1510            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1511               DO jj = 2, jpjm1
1512                  DO ji = 2, jpim1   ! NO vector opt.
1513                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1514                        &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1515                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1516                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1517                        &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1518                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1519                  END DO
1520               END DO
1521            END SELECT
1522         END SELECT
1523         CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
1524         !
1525         !
1526         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
1527            !                                                                     ! Ocean component
1528            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
1529            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
1530            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
1531            zoty1(:,:) = ztmp2(:,:)
1532            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
1533               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
1534               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
1535               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
1536               zity1(:,:) = ztmp2(:,:)
1537            ENDIF
1538         ENDIF
1539         !
1540         ! spherical coordinates to cartesian -> 2 components to 3 components
1541         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
1542            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
1543            ztmp2(:,:) = zoty1(:,:)
1544            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
1545            !
1546            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
1547               ztmp1(:,:) = zitx1(:,:)
1548               ztmp1(:,:) = zity1(:,:)
1549               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
1550            ENDIF
1551         ENDIF
1552         !
1553         IF( ssnd(jps_ocx1)%laction )   CALL cpl_prism_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
1554         IF( ssnd(jps_ocy1)%laction )   CALL cpl_prism_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
1555         IF( ssnd(jps_ocz1)%laction )   CALL cpl_prism_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
1556         !
1557         IF( ssnd(jps_ivx1)%laction )   CALL cpl_prism_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
1558         IF( ssnd(jps_ivy1)%laction )   CALL cpl_prism_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
1559         IF( ssnd(jps_ivz1)%laction )   CALL cpl_prism_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
1560         !
1561      ENDIF
1562      !
1563      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1564      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1565      !
1566      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
1567      !
1568   END SUBROUTINE sbc_cpl_snd
1569   
1570#else
1571   !!----------------------------------------------------------------------
1572   !!   Dummy module                                            NO coupling
1573   !!----------------------------------------------------------------------
1574   USE par_kind        ! kind definition
1575CONTAINS
1576   SUBROUTINE sbc_cpl_snd( kt )
1577      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt
1578   END SUBROUTINE sbc_cpl_snd
1579   !
1580   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
1581      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt, k_fsbc, k_ice
1582   END SUBROUTINE sbc_cpl_rcv
1583   !
1584   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1585      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1586      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1587      p_taui(:,:) = 0.   ;   p_tauj(:,:) = 0. ! stupid definition to avoid warning message when compiling...
1588      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?'
1589   END SUBROUTINE sbc_cpl_ice_tau
1590   !
1591   SUBROUTINE sbc_cpl_ice_flx( p_frld , palbi   , psst    , pist  )
1592      REAL(wp), INTENT(in   ), DIMENSION(:,:  ) ::   p_frld     ! lead fraction                [0 to 1]
1593      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1594      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature      [Celcius]
1595      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature      [Kelvin]
1596      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', p_frld(1,1), palbi(1,1,1), psst(1,1), pist(1,1,1) 
1597   END SUBROUTINE sbc_cpl_ice_flx
1598   
1599#endif
1600
1601   !!======================================================================
1602END MODULE sbccpl