source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 10 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo

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