source: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbccpl.F90 @ 2132

Last change on this file since 2132 was 2132, checked in by omamce, 8 years ago

O.M.

Quick and dirty (very dirty ...) patch to overcome the non allocation of fr*_i0.
See ticket #1139 on NEMO Trac for details.

A better solution is needed in the future.

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