source: branches/UKMO/r5518_sst_landsea_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 7147

Last change on this file since 7147 was 7147, checked in by jcastill, 4 years ago

Remove svn keywords

File size: 120.7 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   !!----------------------------------------------------------------------
12   !!   namsbc_cpl      : coupled formulation namlist
13   !!   sbc_cpl_init    : initialisation of the coupled exchanges
14   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
15   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
16   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
17   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
18   !!   sbc_cpl_snd     : send     fields to the atmosphere
19   !!----------------------------------------------------------------------
20   USE dom_oce         ! ocean space and time domain
21   USE sbc_oce         ! Surface boundary condition: ocean fields
22   USE sbc_ice         ! Surface boundary condition: ice fields
23   USE sbcapr
24   USE sbcdcy          ! surface boundary condition: diurnal cycle
25   USE phycst          ! physical constants
26#if defined key_lim3
27   USE ice             ! ice variables
28#endif
29#if defined key_lim2
30   USE par_ice_2       ! ice parameters
31   USE ice_2           ! ice variables
32#endif
33   USE cpl_oasis3      ! OASIS3 coupling
34   USE geo2ocean       !
35   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
36   USE albedo          !
37   USE in_out_manager  ! I/O manager
38   USE iom             ! NetCDF library
39   USE lib_mpp         ! distribued memory computing library
40   USE wrk_nemo        ! work arrays
41   USE timing          ! Timing
42   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
43   USE eosbn2
44   USE sbcrnf   , ONLY : l_rnfcpl
45#if defined key_cpl_carbon_cycle
46   USE p4zflx, ONLY : oce_co2
47#endif
48#if defined key_cice
49   USE ice_domain_size, only: ncat
50#endif
51#if defined key_lim3
52   USE limthd_dh       ! for CALL lim_thd_snwblow
53#endif
54
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC   sbc_cpl_init       ! routine called by sbcmod.F90
59   PUBLIC   sbc_cpl_rcv        ! routine called by sbc_ice_lim(_2).F90
60   PUBLIC   sbc_cpl_snd        ! routine called by step.F90
61   PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
62   PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).F90
63   PUBLIC   sbc_cpl_alloc      ! routine called in sbcice_cice.F90
64
65   INTEGER, PARAMETER ::   jpr_otx1   =  1            ! 3 atmosphere-ocean stress components on grid 1
66   INTEGER, PARAMETER ::   jpr_oty1   =  2            !
67   INTEGER, PARAMETER ::   jpr_otz1   =  3            !
68   INTEGER, PARAMETER ::   jpr_otx2   =  4            ! 3 atmosphere-ocean stress components on grid 2
69   INTEGER, PARAMETER ::   jpr_oty2   =  5            !
70   INTEGER, PARAMETER ::   jpr_otz2   =  6            !
71   INTEGER, PARAMETER ::   jpr_itx1   =  7            ! 3 atmosphere-ice   stress components on grid 1
72   INTEGER, PARAMETER ::   jpr_ity1   =  8            !
73   INTEGER, PARAMETER ::   jpr_itz1   =  9            !
74   INTEGER, PARAMETER ::   jpr_itx2   = 10            ! 3 atmosphere-ice   stress components on grid 2
75   INTEGER, PARAMETER ::   jpr_ity2   = 11            !
76   INTEGER, PARAMETER ::   jpr_itz2   = 12            !
77   INTEGER, PARAMETER ::   jpr_qsroce = 13            ! Qsr above the ocean
78   INTEGER, PARAMETER ::   jpr_qsrice = 14            ! Qsr above the ice
79   INTEGER, PARAMETER ::   jpr_qsrmix = 15 
80   INTEGER, PARAMETER ::   jpr_qnsoce = 16            ! Qns above the ocean
81   INTEGER, PARAMETER ::   jpr_qnsice = 17            ! Qns above the ice
82   INTEGER, PARAMETER ::   jpr_qnsmix = 18
83   INTEGER, PARAMETER ::   jpr_rain   = 19            ! total liquid precipitation (rain)
84   INTEGER, PARAMETER ::   jpr_snow   = 20            ! solid precipitation over the ocean (snow)
85   INTEGER, PARAMETER ::   jpr_tevp   = 21            ! total evaporation
86   INTEGER, PARAMETER ::   jpr_ievp   = 22            ! solid evaporation (sublimation)
87   INTEGER, PARAMETER ::   jpr_sbpr   = 23            ! sublimation - liquid precipitation - solid precipitation
88   INTEGER, PARAMETER ::   jpr_semp   = 24            ! solid freshwater budget (sublimation - snow)
89   INTEGER, PARAMETER ::   jpr_oemp   = 25            ! ocean freshwater budget (evap - precip)
90   INTEGER, PARAMETER ::   jpr_w10m   = 26            ! 10m wind
91   INTEGER, PARAMETER ::   jpr_dqnsdt = 27            ! d(Q non solar)/d(temperature)
92   INTEGER, PARAMETER ::   jpr_rnf    = 28            ! runoffs
93   INTEGER, PARAMETER ::   jpr_cal    = 29            ! calving
94   INTEGER, PARAMETER ::   jpr_taum   = 30            ! wind stress module
95   INTEGER, PARAMETER ::   jpr_co2    = 31
96   INTEGER, PARAMETER ::   jpr_topm   = 32            ! topmeltn
97   INTEGER, PARAMETER ::   jpr_botm   = 33            ! botmeltn
98   INTEGER, PARAMETER ::   jpr_sflx   = 34            ! salt flux
99   INTEGER, PARAMETER ::   jpr_toce   = 35            ! ocean temperature
100   INTEGER, PARAMETER ::   jpr_soce   = 36            ! ocean salinity
101   INTEGER, PARAMETER ::   jpr_ocx1   = 37            ! ocean current on grid 1
102   INTEGER, PARAMETER ::   jpr_ocy1   = 38            !
103   INTEGER, PARAMETER ::   jpr_ssh    = 39            ! sea surface height
104   INTEGER, PARAMETER ::   jpr_fice   = 40            ! ice fraction         
105   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness
106   INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level
107   INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received
108
109   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
110   INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
111   INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
112   INTEGER, PARAMETER ::   jps_tmix   =  4            ! mixed temperature (ocean+ice)
113   INTEGER, PARAMETER ::   jps_albice =  5            ! ice   albedo
114   INTEGER, PARAMETER ::   jps_albmix =  6            ! mixed albedo
115   INTEGER, PARAMETER ::   jps_hice   =  7            ! ice  thickness
116   INTEGER, PARAMETER ::   jps_hsnw   =  8            ! snow thickness
117   INTEGER, PARAMETER ::   jps_ocx1   =  9            ! ocean current on grid 1
118   INTEGER, PARAMETER ::   jps_ocy1   = 10            !
119   INTEGER, PARAMETER ::   jps_ocz1   = 11            !
120   INTEGER, PARAMETER ::   jps_ivx1   = 12            ! ice   current on grid 1
121   INTEGER, PARAMETER ::   jps_ivy1   = 13            !
122   INTEGER, PARAMETER ::   jps_ivz1   = 14            !
123   INTEGER, PARAMETER ::   jps_co2    = 15
124   INTEGER, PARAMETER ::   jps_soce   = 16            ! ocean salinity
125   INTEGER, PARAMETER ::   jps_ssh    = 17            ! sea surface height
126   INTEGER, PARAMETER ::   jps_qsroce = 18            ! Qsr above the ocean
127   INTEGER, PARAMETER ::   jps_qnsoce = 19            ! Qns above the ocean
128   INTEGER, PARAMETER ::   jps_oemp   = 20            ! ocean freshwater budget (evap - precip)
129   INTEGER, PARAMETER ::   jps_sflx   = 21            ! salt flux
130   INTEGER, PARAMETER ::   jps_otx1   = 22            ! 2 atmosphere-ocean stress components on grid 1
131   INTEGER, PARAMETER ::   jps_oty1   = 23            !
132   INTEGER, PARAMETER ::   jps_rnf    = 24            ! runoffs
133   INTEGER, PARAMETER ::   jps_taum   = 25            ! wind stress module
134   INTEGER, PARAMETER ::   jps_fice2  = 26            ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling)
135   INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl)
136   INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level
137   INTEGER, PARAMETER ::   jpsnd      = 28            ! total number of fields sended
138
139   !                                                         !!** namelist namsbc_cpl **
140   TYPE ::   FLD_C
141      CHARACTER(len = 32) ::   cldes                  ! desciption of the coupling strategy
142      CHARACTER(len = 32) ::   clcat                  ! multiple ice categories strategy
143      CHARACTER(len = 32) ::   clvref                 ! reference of vector ('spherical' or 'cartesian')
144      CHARACTER(len = 32) ::   clvor                  ! orientation of vector fields ('eastward-northward' or 'local grid')
145      CHARACTER(len = 32) ::   clvgrd                 ! grids on which is located the vector fields
146   END TYPE FLD_C
147   ! Send to the atmosphere                           !
148   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                       
149   ! Received from the atmosphere                     !
150   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
151   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                       
152   ! Other namelist parameters                        !
153   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
154   LOGICAL     ::   ln_usecplmask          !  use a coupling mask file to merge data received from several models
155                                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
156   TYPE ::   DYNARR     
157      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
158   END TYPE DYNARR
159
160   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                      ! all fields recieved from the atmosphere
161
162   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
163
164   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
165
166   !! Substitution
167#  include "domzgr_substitute.h90"
168#  include "vectopt_loop_substitute.h90"
169   !!----------------------------------------------------------------------
170   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
171   !! $Id$
172   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
173   !!----------------------------------------------------------------------
174
175CONTAINS
176 
177   INTEGER FUNCTION sbc_cpl_alloc()
178      !!----------------------------------------------------------------------
179      !!             ***  FUNCTION sbc_cpl_alloc  ***
180      !!----------------------------------------------------------------------
181      INTEGER :: ierr(3)
182      !!----------------------------------------------------------------------
183      ierr(:) = 0
184      !
185      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
186     
187#if ! defined key_lim3 && ! defined key_lim2 && ! defined key_cice
188      ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
189#endif
190      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
191      !
192      sbc_cpl_alloc = MAXVAL( ierr )
193      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
194      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
195      !
196   END FUNCTION sbc_cpl_alloc
197
198
199   SUBROUTINE sbc_cpl_init( k_ice )     
200      !!----------------------------------------------------------------------
201      !!             ***  ROUTINE sbc_cpl_init  ***
202      !!
203      !! ** Purpose :   Initialisation of send and received information from
204      !!                the atmospheric component
205      !!
206      !! ** Method  : * Read namsbc_cpl namelist
207      !!              * define the receive interface
208      !!              * define the send    interface
209      !!              * initialise the OASIS coupler
210      !!----------------------------------------------------------------------
211      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
212      !!
213      INTEGER ::   jn   ! dummy loop index
214      INTEGER ::   ios  ! Local integer output status for namelist read
215      INTEGER ::   inum 
216      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
217      !!
218      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      &
219         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      &
220         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   &
221         &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask
222      !!---------------------------------------------------------------------
223      !
224      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
225      !
226      CALL wrk_alloc( jpi,jpj, zacs, zaos )
227
228      ! ================================ !
229      !      Namelist informations       !
230      ! ================================ !
231
232      REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
233      READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
234901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
235
236      REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
237      READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
238902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
239      IF(lwm) WRITE ( numond, namsbc_cpl )
240
241      IF(lwp) THEN                        ! control print
242         WRITE(numout,*)
243         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
244         WRITE(numout,*)'~~~~~~~~~~~~'
245      ENDIF
246      IF( lwp .AND. ln_cpl ) THEN                        ! control print
247         WRITE(numout,*)'  received fields (mutiple ice categogies)'
248         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
249         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
250         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
251         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
252         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
253         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
254         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
255         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
256         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
257         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
258         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
259         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
260         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
261         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
262         WRITE(numout,*)'  sent fields (multiple ice categories)'
263         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
264         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
265         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
266         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
267         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
268         WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
269         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
270         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
271         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
272         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
273      ENDIF
274
275      !                                   ! allocate sbccpl arrays
276      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
277     
278      ! ================================ !
279      !   Define the receive interface   !
280      ! ================================ !
281      nrcvinfo(:) = OASIS_idle   ! needed by nrcvinfo(jpr_otx1) if we do not receive ocean stress
282
283      ! for each field: define the OASIS name                              (srcv(:)%clname)
284      !                 define receive or not from the namelist parameters (srcv(:)%laction)
285      !                 define the north fold type of lbc                  (srcv(:)%nsgn)
286
287      ! default definitions of srcv
288      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
289
290      !                                                      ! ------------------------- !
291      !                                                      ! ice and ocean wind stress !   
292      !                                                      ! ------------------------- !
293      !                                                           ! Name
294      srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
295      srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
296      srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
297      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
298      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
299      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
300      !
301      srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
302      srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
303      srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
304      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
305      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
306      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
307      !
308      ! Vectors: change of sign at north fold ONLY if on the local grid
309      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
310     
311      !                                                           ! Set grid and action
312      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
313      CASE( 'T' ) 
314         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
315         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
316         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
317      CASE( 'U,V' ) 
318         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
319         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
320         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
321         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
322         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
323      CASE( 'U,V,T' )
324         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
325         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
326         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
327         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
328         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
329      CASE( 'U,V,I' )
330         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
331         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
332         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
333         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
334         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
335      CASE( 'U,V,F' )
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  = 'F'        ! ice components given at F-point
339         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
340         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
341      CASE( 'T,I' ) 
342         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
343         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
344         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
345         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
346      CASE( 'T,F' ) 
347         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
348         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
349         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
350         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
351      CASE( 'T,U,V' )
352         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
353         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
354         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
355         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
356         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
357      CASE default   
358         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
359      END SELECT
360      !
361      IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
362         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE. 
363      !
364      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
365            srcv(jpr_otx2:jpr_otz2)%laction = .FALSE. 
366            srcv(jpr_itx2:jpr_itz2)%laction = .FALSE. 
367            srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid   ! not needed but cleaner...
368            srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid   ! not needed but cleaner...
369      ENDIF
370      !
371      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
372         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
373         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
374         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
375      ENDIF
376       
377      !                                                      ! ------------------------- !
378      !                                                      !    freshwater budget      !   E-P
379      !                                                      ! ------------------------- !
380      ! we suppose that atmosphere modele do not make the difference between precipiration (liquide or solid)
381      ! over ice of free ocean within the same atmospheric cell.cd
382      srcv(jpr_rain)%clname = 'OTotRain'      ! Rain = liquid precipitation
383      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
384      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
385      srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation
386      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation
387      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
388      srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
389      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
390      CASE( 'none'          )       ! nothing to do
391      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
392      CASE( 'conservative'  )
393         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
394         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
395      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
396      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
397      END SELECT
398
399      !                                                      ! ------------------------- !
400      !                                                      !     Runoffs & Calving     !   
401      !                                                      ! ------------------------- !
402      srcv(jpr_rnf   )%clname = 'O_Runoff'
403      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
404         srcv(jpr_rnf)%laction = .TRUE.
405         l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
406         ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
407         IF(lwp) WRITE(numout,*)
408         IF(lwp) WRITE(numout,*) '   runoffs received from oasis -> force ln_rnf = ', ln_rnf
409      ENDIF
410      !
411      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
412
413      !                                                      ! ------------------------- !
414      !                                                      !    non solar radiation    !   Qns
415      !                                                      ! ------------------------- !
416      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
417      srcv(jpr_qnsice)%clname = 'O_QnsIce'
418      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
419      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
420      CASE( 'none'          )       ! nothing to do
421      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
422      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
423      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
424      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
425      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
426      END SELECT
427      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
428         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
429      !                                                      ! ------------------------- !
430      !                                                      !    solar radiation        !   Qsr
431      !                                                      ! ------------------------- !
432      srcv(jpr_qsroce)%clname = 'O_QsrOce'
433      srcv(jpr_qsrice)%clname = 'O_QsrIce'
434      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
435      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
436      CASE( 'none'          )       ! nothing to do
437      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
438      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
439      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
440      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
441      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
442      END SELECT
443      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
444         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
445      !                                                      ! ------------------------- !
446      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
447      !                                                      ! ------------------------- !
448      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
449      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
450      !
451      ! non solar sensitivity mandatory for LIM ice model
452      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4 .AND. nn_components /= jp_iam_sas ) &
453         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
454      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
455      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
456         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
457      !                                                      ! ------------------------- !
458      !                                                      !      10m wind module      !   
459      !                                                      ! ------------------------- !
460      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
461      !
462      !                                                      ! ------------------------- !
463      !                                                      !   wind stress module      !   
464      !                                                      ! ------------------------- !
465      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
466      lhftau = srcv(jpr_taum)%laction
467
468      !                                                      ! ------------------------- !
469      !                                                      !      Atmospheric CO2      !
470      !                                                      ! ------------------------- !
471      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
472      !                                                      ! ------------------------- !
473      !                                                      !   topmelt and botmelt     !   
474      !                                                      ! ------------------------- !
475      srcv(jpr_topm )%clname = 'OTopMlt'
476      srcv(jpr_botm )%clname = 'OBotMlt'
477      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
478         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
479            srcv(jpr_topm:jpr_botm)%nct = jpl
480         ELSE
481            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
482         ENDIF
483         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
484      ENDIF
485      !                                                      ! ------------------------------- !
486      !                                                      !   OPA-SAS coupling - rcv by opa !   
487      !                                                      ! ------------------------------- !
488      srcv(jpr_sflx)%clname = 'O_SFLX'
489      srcv(jpr_fice)%clname = 'RIceFrc'
490      !
491      IF( nn_components == jp_iam_opa ) THEN    ! OPA coupled to SAS via OASIS: force received field by OPA (sent by SAS)
492         srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
493         srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
494         srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
495         srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
496         srcv(jpr_otx1)%clgrid = 'U'        ! oce components given at U-point
497         srcv(jpr_oty1)%clgrid = 'V'        !           and           V-point
498         ! Vectors: change of sign at north fold ONLY if on the local grid
499         srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
500         sn_rcv_tau%clvgrd = 'U,V'
501         sn_rcv_tau%clvor = 'local grid'
502         sn_rcv_tau%clvref = 'spherical'
503         sn_rcv_emp%cldes = 'oce only'
504         !
505         IF(lwp) THEN                        ! control print
506            WRITE(numout,*)
507            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
508            WRITE(numout,*)'               OPA component  '
509            WRITE(numout,*)
510            WRITE(numout,*)'  received fields from SAS component '
511            WRITE(numout,*)'                  ice cover '
512            WRITE(numout,*)'                  oce only EMP  '
513            WRITE(numout,*)'                  salt flux  '
514            WRITE(numout,*)'                  mixed oce-ice solar flux  '
515            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
516            WRITE(numout,*)'                  wind stress U,V on local grid and sperical coordinates '
517            WRITE(numout,*)'                  wind stress module'
518            WRITE(numout,*)
519         ENDIF
520      ENDIF
521      !                                                      ! -------------------------------- !
522      !                                                      !   OPA-SAS coupling - rcv by sas  !   
523      !                                                      ! -------------------------------- !
524      srcv(jpr_toce  )%clname = 'I_SSTSST'
525      srcv(jpr_soce  )%clname = 'I_SSSal'
526      srcv(jpr_ocx1  )%clname = 'I_OCurx1'
527      srcv(jpr_ocy1  )%clname = 'I_OCury1'
528      srcv(jpr_ssh   )%clname = 'I_SSHght'
529      srcv(jpr_e3t1st)%clname = 'I_E3T1st'   
530      srcv(jpr_fraqsr)%clname = 'I_FraQsr'   
531      !
532      IF( nn_components == jp_iam_sas ) THEN
533         IF( .NOT. ln_cpl ) srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
534         IF( .NOT. ln_cpl ) srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
535         IF( .NOT. ln_cpl ) srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
536         srcv( (/jpr_toce, jpr_soce, jpr_ssh, jpr_fraqsr, jpr_ocx1, jpr_ocy1/) )%laction = .TRUE.
537         srcv( jpr_e3t1st )%laction = lk_vvl
538         srcv(jpr_ocx1)%clgrid = 'U'        ! oce components given at U-point
539         srcv(jpr_ocy1)%clgrid = 'V'        !           and           V-point
540         ! Vectors: change of sign at north fold ONLY if on the local grid
541         srcv(jpr_ocx1:jpr_ocy1)%nsgn = -1.
542         ! Change first letter to couple with atmosphere if already coupled OPA
543         ! this is nedeed as each variable name used in the namcouple must be unique:
544         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
545         DO jn = 1, jprcv
546            IF ( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
547         END DO
548         !
549         IF(lwp) THEN                        ! control print
550            WRITE(numout,*)
551            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
552            WRITE(numout,*)'               SAS component  '
553            WRITE(numout,*)
554            IF( .NOT. ln_cpl ) THEN
555               WRITE(numout,*)'  received fields from OPA component '
556            ELSE
557               WRITE(numout,*)'  Additional received fields from OPA component : '
558            ENDIF
559            WRITE(numout,*)'               sea surface temperature (Celcius) '
560            WRITE(numout,*)'               sea surface salinity ' 
561            WRITE(numout,*)'               surface currents ' 
562            WRITE(numout,*)'               sea surface height ' 
563            WRITE(numout,*)'               thickness of first ocean T level '       
564            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
565            WRITE(numout,*)
566         ENDIF
567      ENDIF
568     
569      ! =================================================== !
570      ! Allocate all parts of frcv used for received fields !
571      ! =================================================== !
572      DO jn = 1, jprcv
573         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
574      END DO
575      ! Allocate taum part of frcv which is used even when not received as coupling field
576      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
577      ! Allocate w10m part of frcv which is used even when not received as coupling field
578      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
579      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
580      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
581      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
582      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
583      IF( k_ice /= 0 ) THEN
584         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
585         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
586      END IF
587
588      ! ================================ !
589      !     Define the send interface    !
590      ! ================================ !
591      ! for each field: define the OASIS name                           (ssnd(:)%clname)
592      !                 define send or not from the namelist parameters (ssnd(:)%laction)
593      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
594     
595      ! default definitions of nsnd
596      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
597         
598      !                                                      ! ------------------------- !
599      !                                                      !    Surface temperature    !
600      !                                                      ! ------------------------- !
601      ssnd(jps_toce)%clname = 'O_SSTSST'
602      ssnd(jps_tice)%clname = 'O_TepIce'
603      ssnd(jps_tmix)%clname = 'O_TepMix'
604      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
605      CASE( 'none'                                 )       ! nothing to do
606      CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE.
607      CASE( 'oce and ice' , 'weighted oce and ice' )
608         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
609         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
610      CASE( 'mixed oce-ice'                        )   ;   ssnd( jps_tmix )%laction = .TRUE.
611      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
612      END SELECT
613           
614      !                                                      ! ------------------------- !
615      !                                                      !          Albedo           !
616      !                                                      ! ------------------------- !
617      ssnd(jps_albice)%clname = 'O_AlbIce' 
618      ssnd(jps_albmix)%clname = 'O_AlbMix'
619      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
620      CASE( 'none'                 )     ! nothing to do
621      CASE( 'ice' , 'weighted ice' )   ; ssnd(jps_albice)%laction = .TRUE.
622      CASE( 'mixed oce-ice'        )   ; ssnd(jps_albmix)%laction = .TRUE.
623      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
624      END SELECT
625      !
626      ! Need to calculate oceanic albedo if
627      !     1. sending mixed oce-ice albedo or
628      !     2. receiving mixed oce-ice solar radiation
629      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
630         CALL albedo_oce( zaos, zacs )
631         ! Due to lack of information on nebulosity : mean clear/overcast sky
632         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
633      ENDIF
634
635      !                                                      ! ------------------------- !
636      !                                                      !  Ice fraction & Thickness !
637      !                                                      ! ------------------------- !
638      ssnd(jps_fice)%clname = 'OIceFrc'
639      ssnd(jps_hice)%clname = 'OIceTck'
640      ssnd(jps_hsnw)%clname = 'OSnwTck'
641      IF( k_ice /= 0 ) THEN
642         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
643! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
644         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
645      ENDIF
646     
647      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
648      CASE( 'none'         )       ! nothing to do
649      CASE( 'ice and snow' ) 
650         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
651         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
652            ssnd(jps_hice:jps_hsnw)%nct = jpl
653         ENDIF
654      CASE ( 'weighted ice and snow' ) 
655         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
656         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
657      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
658      END SELECT
659
660      !                                                      ! ------------------------- !
661      !                                                      !      Surface current      !
662      !                                                      ! ------------------------- !
663      !        ocean currents              !            ice velocities
664      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
665      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
666      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
667      !
668      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
669
670      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
671         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
672      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
673         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
674         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
675      ENDIF
676      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
677      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
678      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
679      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
680      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
681      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
682      CASE( 'weighted oce and ice' )   !   nothing to do
683      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
684      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
685      END SELECT
686
687      !                                                      ! ------------------------- !
688      !                                                      !          CO2 flux         !
689      !                                                      ! ------------------------- !
690      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
691
692      !                                                      ! ------------------------------- !
693      !                                                      !   OPA-SAS coupling - snd by opa !   
694      !                                                      ! ------------------------------- !
695      ssnd(jps_ssh   )%clname = 'O_SSHght' 
696      ssnd(jps_soce  )%clname = 'O_SSSal' 
697      ssnd(jps_e3t1st)%clname = 'O_E3T1st'   
698      ssnd(jps_fraqsr)%clname = 'O_FraQsr'
699      !
700      IF( nn_components == jp_iam_opa ) THEN
701         ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
702         ssnd( (/jps_toce, jps_soce, jps_ssh, jps_fraqsr, jps_ocx1, jps_ocy1/) )%laction = .TRUE.
703         ssnd( jps_e3t1st )%laction = lk_vvl
704         ! vector definition: not used but cleaner...
705         ssnd(jps_ocx1)%clgrid  = 'U'        ! oce components given at U-point
706         ssnd(jps_ocy1)%clgrid  = 'V'        !           and           V-point
707         sn_snd_crt%clvgrd = 'U,V'
708         sn_snd_crt%clvor = 'local grid'
709         sn_snd_crt%clvref = 'spherical'
710         !
711         IF(lwp) THEN                        ! control print
712            WRITE(numout,*)
713            WRITE(numout,*)'  sent fields to SAS component '
714            WRITE(numout,*)'               sea surface temperature (T before, Celcius) '
715            WRITE(numout,*)'               sea surface salinity ' 
716            WRITE(numout,*)'               surface currents U,V on local grid and spherical coordinates' 
717            WRITE(numout,*)'               sea surface height ' 
718            WRITE(numout,*)'               thickness of first ocean T level '       
719            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
720            WRITE(numout,*)
721         ENDIF
722      ENDIF
723      !                                                      ! ------------------------------- !
724      !                                                      !   OPA-SAS coupling - snd by sas !   
725      !                                                      ! ------------------------------- !
726      ssnd(jps_sflx  )%clname = 'I_SFLX'     
727      ssnd(jps_fice2 )%clname = 'IIceFrc'
728      ssnd(jps_qsroce)%clname = 'I_QsrOce'   
729      ssnd(jps_qnsoce)%clname = 'I_QnsOce'   
730      ssnd(jps_oemp  )%clname = 'IOEvaMPr' 
731      ssnd(jps_otx1  )%clname = 'I_OTaux1'   
732      ssnd(jps_oty1  )%clname = 'I_OTauy1'   
733      ssnd(jps_rnf   )%clname = 'I_Runoff'   
734      ssnd(jps_taum  )%clname = 'I_TauMod'   
735      !
736      IF( nn_components == jp_iam_sas ) THEN
737         IF( .NOT. ln_cpl ) ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
738         ssnd( (/jps_qsroce, jps_qnsoce, jps_oemp, jps_fice2, jps_sflx, jps_otx1, jps_oty1, jps_taum/) )%laction = .TRUE.
739         !
740         ! Change first letter to couple with atmosphere if already coupled with sea_ice
741         ! this is nedeed as each variable name used in the namcouple must be unique:
742         ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere
743         DO jn = 1, jpsnd
744            IF ( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname))
745         END DO
746         !
747         IF(lwp) THEN                        ! control print
748            WRITE(numout,*)
749            IF( .NOT. ln_cpl ) THEN
750               WRITE(numout,*)'  sent fields to OPA component '
751            ELSE
752               WRITE(numout,*)'  Additional sent fields to OPA component : '
753            ENDIF
754            WRITE(numout,*)'                  ice cover '
755            WRITE(numout,*)'                  oce only EMP  '
756            WRITE(numout,*)'                  salt flux  '
757            WRITE(numout,*)'                  mixed oce-ice solar flux  '
758            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
759            WRITE(numout,*)'                  wind stress U,V components'
760            WRITE(numout,*)'                  wind stress module'
761         ENDIF
762      ENDIF
763
764      !
765      ! ================================ !
766      !   initialisation of the coupler  !
767      ! ================================ !
768
769      CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
770     
771      IF (ln_usecplmask) THEN
772         xcplmask(:,:,:) = 0.
773         CALL iom_open( 'cplmask', inum )
774         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:nlci,1:nlcj,1:nn_cplmodel),   &
775            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nn_cplmodel /) )
776         CALL iom_close( inum )
777      ELSE
778         xcplmask(:,:,:) = 1.
779      ENDIF
780      xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
781      !
782      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
783      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
784         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
785      ncpl_qsr_freq = 86400 / ncpl_qsr_freq
786
787      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
788      !
789      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
790      !
791   END SUBROUTINE sbc_cpl_init
792
793
794   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
795      !!----------------------------------------------------------------------
796      !!             ***  ROUTINE sbc_cpl_rcv  ***
797      !!
798      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
799      !!                provide the ocean heat and freshwater fluxes.
800      !!
801      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
802      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
803      !!                to know if the field was really received or not
804      !!
805      !!              --> If ocean stress was really received:
806      !!
807      !!                  - transform the received ocean stress vector from the received
808      !!                 referential and grid into an atmosphere-ocean stress in
809      !!                 the (i,j) ocean referencial and at the ocean velocity point.
810      !!                    The received stress are :
811      !!                     - defined by 3 components (if cartesian coordinate)
812      !!                            or by 2 components (if spherical)
813      !!                     - oriented along geographical   coordinate (if eastward-northward)
814      !!                            or  along the local grid coordinate (if local grid)
815      !!                     - given at U- and V-point, resp.   if received on 2 grids
816      !!                            or at T-point               if received on 1 grid
817      !!                    Therefore and if necessary, they are successively
818      !!                  processed in order to obtain them
819      !!                     first  as  2 components on the sphere
820      !!                     second as  2 components oriented along the local grid
821      !!                     third  as  2 components on the U,V grid
822      !!
823      !!              -->
824      !!
825      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
826      !!             and total ocean freshwater fluxes 
827      !!
828      !! ** Method  :   receive all fields from the atmosphere and transform
829      !!              them into ocean surface boundary condition fields
830      !!
831      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
832      !!                        taum         wind stress module at T-point
833      !!                        wndm         wind speed  module at T-point over free ocean or leads in presence of sea-ice
834      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
835      !!                                     and the latent heat flux of solid precip. melting
836      !!                        qsr          solar ocean heat fluxes   (ocean only case)
837      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
838      !!----------------------------------------------------------------------
839      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
840      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation
841      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
842
843      !!
844      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
845      INTEGER  ::   ji, jj, jn             ! dummy loop indices
846      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
847      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
848      REAL(wp) ::   zcoef                  ! temporary scalar
849      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
850      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
851      REAL(wp) ::   zzx, zzy               ! temporary variables
852      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr
853      !!----------------------------------------------------------------------
854      !
855      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
856      !
857      CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
858      !
859      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
860      !
861      !                                                      ! ======================================================= !
862      !                                                      ! Receive all the atmos. fields (including ice information)
863      !                                                      ! ======================================================= !
864      isec = ( kt - nit000 ) * NINT( rdttra(1) )                ! date of exchanges
865      DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
866         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
867      END DO
868
869      !                                                      ! ========================= !
870      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
871         !                                                   ! ========================= !
872         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
873         ! => need to be done only when we receive the field
874         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
875            !
876            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
877               !                                                       ! (cartesian to spherical -> 3 to 2 components)
878               !
879               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
880                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
881               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
882               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
883               !
884               IF( srcv(jpr_otx2)%laction ) THEN
885                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
886                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
887                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
888                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
889               ENDIF
890               !
891            ENDIF
892            !
893            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
894               !                                                       ! (geographical to local grid -> rotate the components)
895               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
896               IF( srcv(jpr_otx2)%laction ) THEN
897                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
898               ELSE 
899                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
900               ENDIF
901               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
902               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
903            ENDIF
904            !                             
905            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
906               DO jj = 2, jpjm1                                          ! T ==> (U,V)
907                  DO ji = fs_2, fs_jpim1   ! vector opt.
908                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
909                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
910                  END DO
911               END DO
912               CALL lbc_lnk( frcv(jpr_otx1)%z3(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
913            ENDIF
914            llnewtx = .TRUE.
915         ELSE
916            llnewtx = .FALSE.
917         ENDIF
918         !                                                   ! ========================= !
919      ELSE                                                   !   No dynamical coupling   !
920         !                                                   ! ========================= !
921         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
922         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
923         llnewtx = .TRUE.
924         !
925      ENDIF
926      !                                                      ! ========================= !
927      !                                                      !    wind stress module     !   (taum)
928      !                                                      ! ========================= !
929      !
930      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
931         ! => need to be done only when otx1 was changed
932         IF( llnewtx ) THEN
933            DO jj = 2, jpjm1
934               DO ji = fs_2, fs_jpim1   ! vect. opt.
935                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
936                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
937                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
938               END DO
939            END DO
940            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
941            llnewtau = .TRUE.
942         ELSE
943            llnewtau = .FALSE.
944         ENDIF
945      ELSE
946         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
947         ! Stress module can be negative when received (interpolation problem)
948         IF( llnewtau ) THEN
949            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
950         ENDIF
951      ENDIF
952      !
953      !                                                      ! ========================= !
954      !                                                      !      10 m wind speed      !   (wndm)
955      !                                                      ! ========================= !
956      !
957      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
958         ! => need to be done only when taumod was changed
959         IF( llnewtau ) THEN
960            zcoef = 1. / ( zrhoa * zcdrag ) 
961            DO jj = 1, jpj
962               DO ji = 1, jpi 
963                  frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
964               END DO
965            END DO
966         ENDIF
967      ENDIF
968
969      ! u(v)tau and taum will be modified by ice model
970      ! -> need to be reset before each call of the ice/fsbc     
971      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
972         !
973         IF( ln_mixcpl ) THEN
974            utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)
975            vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
976            taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
977            wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
978         ELSE
979            utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
980            vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
981            taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
982            wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
983         ENDIF
984         CALL iom_put( "taum_oce", taum )   ! output wind stress module
985         
986      ENDIF
987
988#if defined key_cpl_carbon_cycle
989      !                                                      ! ================== !
990      !                                                      ! atmosph. CO2 (ppm) !
991      !                                                      ! ================== !
992      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
993#endif
994
995      !  Fields received by SAS when OASIS coupling
996      !  (arrays no more filled at sbcssm stage)
997      !                                                      ! ================== !
998      !                                                      !        SSS         !
999      !                                                      ! ================== !
1000      IF( srcv(jpr_soce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1001         sss_m(:,:) = frcv(jpr_soce)%z3(:,:,1)
1002         CALL iom_put( 'sss_m', sss_m )
1003      ENDIF
1004      !                                               
1005      !                                                      ! ================== !
1006      !                                                      !        SST         !
1007      !                                                      ! ================== !
1008      IF( srcv(jpr_toce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1009         sst_m(:,:) = frcv(jpr_toce)%z3(:,:,1)
1010         IF( srcv(jpr_soce)%laction .AND. ln_useCT ) THEN    ! make sure that sst_m is the potential temperature
1011            sst_m(:,:) = eos_pt_from_ct( sst_m(:,:), sss_m(:,:) )
1012         ENDIF
1013      ENDIF
1014      !                                                      ! ================== !
1015      !                                                      !        SSH         !
1016      !                                                      ! ================== !
1017      IF( srcv(jpr_ssh )%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1018         ssh_m(:,:) = frcv(jpr_ssh )%z3(:,:,1)
1019         CALL iom_put( 'ssh_m', ssh_m )
1020      ENDIF
1021      !                                                      ! ================== !
1022      !                                                      !  surface currents  !
1023      !                                                      ! ================== !
1024      IF( srcv(jpr_ocx1)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1025         ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
1026         ub (:,:,1) = ssu_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
1027         CALL iom_put( 'ssu_m', ssu_m )
1028      ENDIF
1029      IF( srcv(jpr_ocy1)%laction ) THEN
1030         ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
1031         vb (:,:,1) = ssv_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
1032         CALL iom_put( 'ssv_m', ssv_m )
1033      ENDIF
1034      !                                                      ! ======================== !
1035      !                                                      !  first T level thickness !
1036      !                                                      ! ======================== !
1037      IF( srcv(jpr_e3t1st )%laction ) THEN                   ! received by sas in case of opa <-> sas coupling
1038         e3t_m(:,:) = frcv(jpr_e3t1st )%z3(:,:,1)
1039         CALL iom_put( 'e3t_m', e3t_m(:,:) )
1040      ENDIF
1041      !                                                      ! ================================ !
1042      !                                                      !  fraction of solar net radiation !
1043      !                                                      ! ================================ !
1044      IF( srcv(jpr_fraqsr)%laction ) THEN                    ! received by sas in case of opa <-> sas coupling
1045         frq_m(:,:) = frcv(jpr_fraqsr)%z3(:,:,1)
1046         CALL iom_put( 'frq_m', frq_m )
1047      ENDIF
1048     
1049      !                                                      ! ========================= !
1050      IF( k_ice <= 1 .AND. MOD( kt-1, k_fsbc ) == 0 ) THEN   !  heat & freshwater fluxes ! (Ocean only case)
1051         !                                                   ! ========================= !
1052         !
1053         !                                                       ! total freshwater fluxes over the ocean (emp)
1054         IF( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction ) THEN
1055            SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
1056            CASE( 'conservative' )
1057               zemp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
1058            CASE( 'oce only', 'oce and ice' )
1059               zemp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
1060            CASE default
1061               CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
1062            END SELECT
1063         ELSE
1064            zemp(:,:) = 0._wp
1065         ENDIF
1066         !
1067         !                                                        ! runoffs and calving (added in emp)
1068         IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1069         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
1070         
1071         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
1072         ELSE                   ;   emp(:,:) =                              zemp(:,:)
1073         ENDIF
1074         !
1075         !                                                       ! non solar heat flux over the ocean (qns)
1076         IF(      srcv(jpr_qnsoce)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
1077         ELSE IF( srcv(jpr_qnsmix)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
1078         ELSE                                       ;   zqns(:,:) = 0._wp
1079         END IF
1080         ! update qns over the free ocean with:
1081         IF( nn_components /= jp_iam_opa ) THEN
1082            zqns(:,:) =  zqns(:,:) - zemp(:,:) * sst_m(:,:) * rcp         ! remove heat content due to mass flux (assumed to be at SST)
1083            IF( srcv(jpr_snow  )%laction ) THEN
1084               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean
1085            ENDIF
1086         ENDIF
1087         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
1088         ELSE                   ;   qns(:,:) =                              zqns(:,:)
1089         ENDIF
1090
1091         !                                                       ! solar flux over the ocean          (qsr)
1092         IF     ( srcv(jpr_qsroce)%laction ) THEN   ;   zqsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
1093         ELSE IF( srcv(jpr_qsrmix)%laction ) then   ;   zqsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
1094         ELSE                                       ;   zqsr(:,:) = 0._wp
1095         ENDIF
1096         IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle
1097         IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)
1098         ELSE                   ;   qsr(:,:) =                              zqsr(:,:)
1099         ENDIF
1100         !
1101         ! salt flux over the ocean (received by opa in case of opa <-> sas coupling)
1102         IF( srcv(jpr_sflx )%laction )   sfx(:,:) = frcv(jpr_sflx  )%z3(:,:,1)
1103         ! Ice cover  (received by opa in case of opa <-> sas coupling)
1104         IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
1105         !
1106
1107      ENDIF
1108      !
1109      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
1110      !
1111      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
1112      !
1113   END SUBROUTINE sbc_cpl_rcv
1114   
1115
1116   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1117      !!----------------------------------------------------------------------
1118      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
1119      !!
1120      !! ** Purpose :   provide the stress over sea-ice in coupled mode
1121      !!
1122      !! ** Method  :   transform the received stress from the atmosphere into
1123      !!             an atmosphere-ice stress in the (i,j) ocean referencial
1124      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
1125      !!                'C'-grid : i- (j-) components given at U- (V-) point
1126      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
1127      !!
1128      !!                The received stress are :
1129      !!                 - defined by 3 components (if cartesian coordinate)
1130      !!                        or by 2 components (if spherical)
1131      !!                 - oriented along geographical   coordinate (if eastward-northward)
1132      !!                        or  along the local grid coordinate (if local grid)
1133      !!                 - given at U- and V-point, resp.   if received on 2 grids
1134      !!                        or at a same point (T or I) if received on 1 grid
1135      !!                Therefore and if necessary, they are successively
1136      !!             processed in order to obtain them
1137      !!                 first  as  2 components on the sphere
1138      !!                 second as  2 components oriented along the local grid
1139      !!                 third  as  2 components on the cp_ice_msh point
1140      !!
1141      !!                Except in 'oce and ice' case, only one vector stress field
1142      !!             is received. It has already been processed in sbc_cpl_rcv
1143      !!             so that it is now defined as (i,j) components given at U-
1144      !!             and V-points, respectively. Therefore, only the third
1145      !!             transformation is done and only if the ice-grid is a 'I'-grid.
1146      !!
1147      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
1148      !!----------------------------------------------------------------------
1149      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1150      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1151      !!
1152      INTEGER ::   ji, jj                          ! dummy loop indices
1153      INTEGER ::   itx                             ! index of taux over ice
1154      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
1155      !!----------------------------------------------------------------------
1156      !
1157      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
1158      !
1159      CALL wrk_alloc( jpi,jpj, ztx, zty )
1160
1161      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
1162      ELSE                                ;   itx =  jpr_otx1
1163      ENDIF
1164
1165      ! do something only if we just received the stress from atmosphere
1166      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
1167
1168         !                                                      ! ======================= !
1169         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
1170            !                                                   ! ======================= !
1171           
1172            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
1173               !                                                       ! (cartesian to spherical -> 3 to 2 components)
1174               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
1175                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
1176               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
1177               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
1178               !
1179               IF( srcv(jpr_itx2)%laction ) THEN
1180                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
1181                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
1182                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
1183                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
1184               ENDIF
1185               !
1186            ENDIF
1187            !
1188            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
1189               !                                                       ! (geographical to local grid -> rotate the components)
1190               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
1191               IF( srcv(jpr_itx2)%laction ) THEN
1192                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
1193               ELSE
1194                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
1195               ENDIF
1196               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1197               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
1198            ENDIF
1199            !                                                   ! ======================= !
1200         ELSE                                                   !     use ocean stress    !
1201            !                                                   ! ======================= !
1202            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
1203            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
1204            !
1205         ENDIF
1206         !                                                      ! ======================= !
1207         !                                                      !     put on ice grid     !
1208         !                                                      ! ======================= !
1209         !   
1210         !                                                  j+1   j     -----V---F
1211         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
1212         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
1213         !                                                               |       |
1214         !                                                   j    j-1   -I-------|
1215         !                                               (for I)         |       |
1216         !                                                              i-1  i   i
1217         !                                                               i      i+1 (for I)
1218         SELECT CASE ( cp_ice_msh )
1219            !
1220         CASE( 'I' )                                         ! B-grid ==> I
1221            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1222            CASE( 'U' )
1223               DO jj = 2, jpjm1                                   ! (U,V) ==> I
1224                  DO ji = 2, jpim1   ! NO vector opt.
1225                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
1226                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1227                  END DO
1228               END DO
1229            CASE( 'F' )
1230               DO jj = 2, jpjm1                                   ! F ==> I
1231                  DO ji = 2, jpim1   ! NO vector opt.
1232                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
1233                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
1234                  END DO
1235               END DO
1236            CASE( 'T' )
1237               DO jj = 2, jpjm1                                   ! T ==> I
1238                  DO ji = 2, jpim1   ! NO vector opt.
1239                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
1240                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
1241                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
1242                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1243                  END DO
1244               END DO
1245            CASE( 'I' )
1246               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
1247               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1248            END SELECT
1249            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1250               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1251            ENDIF
1252            !
1253         CASE( 'F' )                                         ! B-grid ==> F
1254            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1255            CASE( 'U' )
1256               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1257                  DO ji = fs_2, fs_jpim1   ! vector opt.
1258                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1259                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1260                  END DO
1261               END DO
1262            CASE( 'I' )
1263               DO jj = 2, jpjm1                                   ! I ==> F
1264                  DO ji = 2, jpim1   ! NO vector opt.
1265                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1266                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1267                  END DO
1268               END DO
1269            CASE( 'T' )
1270               DO jj = 2, jpjm1                                   ! T ==> F
1271                  DO ji = 2, jpim1   ! NO vector opt.
1272                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1273                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1274                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1275                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1276                  END DO
1277               END DO
1278            CASE( 'F' )
1279               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1280               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1281            END SELECT
1282            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1283               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1284            ENDIF
1285            !
1286         CASE( 'C' )                                         ! C-grid ==> U,V
1287            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1288            CASE( 'U' )
1289               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1290               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1291            CASE( 'F' )
1292               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1293                  DO ji = fs_2, fs_jpim1   ! vector opt.
1294                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1295                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1296                  END DO
1297               END DO
1298            CASE( 'T' )
1299               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1300                  DO ji = fs_2, fs_jpim1   ! vector opt.
1301                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1302                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1303                  END DO
1304               END DO
1305            CASE( 'I' )
1306               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1307                  DO ji = 2, jpim1   ! NO vector opt.
1308                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1309                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1310                  END DO
1311               END DO
1312            END SELECT
1313            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1314               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1315            ENDIF
1316         END SELECT
1317
1318      ENDIF
1319      !   
1320      CALL wrk_dealloc( jpi,jpj, ztx, zty )
1321      !
1322      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
1323      !
1324   END SUBROUTINE sbc_cpl_ice_tau
1325   
1326
1327   SUBROUTINE sbc_cpl_ice_flx( p_frld, palbi, psst, pist )
1328      !!----------------------------------------------------------------------
1329      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1330      !!
1331      !! ** Purpose :   provide the heat and freshwater fluxes of the
1332      !!              ocean-ice system.
1333      !!
1334      !! ** Method  :   transform the fields received from the atmosphere into
1335      !!             surface heat and fresh water boundary condition for the
1336      !!             ice-ocean system. The following fields are provided:
1337      !!              * total non solar, solar and freshwater fluxes (qns_tot,
1338      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1339      !!             NB: emp_tot include runoffs and calving.
1340      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1341      !!             emp_ice = sublimation - solid precipitation as liquid
1342      !!             precipitation are re-routed directly to the ocean and
1343      !!             runoffs and calving directly enter the ocean.
1344      !!              * solid precipitation (sprecip), used to add to qns_tot
1345      !!             the heat lost associated to melting solid precipitation
1346      !!             over the ocean fraction.
1347      !!       ===>> CAUTION here this changes the net heat flux received from
1348      !!             the atmosphere
1349      !!
1350      !!                  - the fluxes have been separated from the stress as
1351      !!                 (a) they are updated at each ice time step compare to
1352      !!                 an update at each coupled time step for the stress, and
1353      !!                 (b) the conservative computation of the fluxes over the
1354      !!                 sea-ice area requires the knowledge of the ice fraction
1355      !!                 after the ice advection and before the ice thermodynamics,
1356      !!                 so that the stress is updated before the ice dynamics
1357      !!                 while the fluxes are updated after it.
1358      !!
1359      !! ** Action  :   update at each nf_ice time step:
1360      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1361      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1362      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
1363      !!                   emp_ice            ice sublimation - solid precipitation over the ice
1364      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
1365      !!                   sprecip             solid precipitation over the ocean 
1366      !!----------------------------------------------------------------------
1367      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1368      ! optional arguments, used only in 'mixed oce-ice' case
1369      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo
1370      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
1371      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
1372      !
1373      INTEGER ::   jl         ! dummy loop index
1374      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk
1375      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
1376      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice
1377      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3
1378      !!----------------------------------------------------------------------
1379      !
1380      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1381      !
1382      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
1383      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
1384
1385      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
1386      zicefr(:,:) = 1.- p_frld(:,:)
1387      zcptn(:,:) = rcp * sst_m(:,:)
1388      !
1389      !                                                      ! ========================= !
1390      !                                                      !    freshwater budget      !   (emp)
1391      !                                                      ! ========================= !
1392      !
1393      !                                                           ! total Precipitation - total Evaporation (emp_tot)
1394      !                                                           ! solid precipitation - sublimation       (emp_ice)
1395      !                                                           ! solid Precipitation                     (sprecip)
1396      !                                                           ! liquid + solid Precipitation            (tprecip)
1397      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1398      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1399         zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
1400         ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
1401         zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
1402         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1403            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation
1404         IF( iom_use('hflx_rain_cea') )   &
1405            CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
1406         IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   &
1407            ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1408         IF( iom_use('evap_ao_cea'  ) )   &
1409            CALL iom_put( 'evap_ao_cea'  , ztmp                   )   ! ice-free oce evap (cell average)
1410         IF( iom_use('hflx_evap_cea') )   &
1411            CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) )   ! heat flux from from evap (cell average)
1412      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1413         zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1414         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
1415         zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
1416         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
1417      END SELECT
1418
1419      IF( iom_use('subl_ai_cea') )   &
1420         CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
1421      !   
1422      !                                                           ! runoffs and calving (put in emp_tot)
1423      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1424      IF( srcv(jpr_cal)%laction ) THEN
1425         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1426         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
1427      ENDIF
1428
1429      IF( ln_mixcpl ) THEN
1430         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
1431         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
1432         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
1433         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
1434      ELSE
1435         emp_tot(:,:) =                                  zemp_tot(:,:)
1436         emp_ice(:,:) =                                  zemp_ice(:,:)
1437         sprecip(:,:) =                                  zsprecip(:,:)
1438         tprecip(:,:) =                                  ztprecip(:,:)
1439      ENDIF
1440
1441         CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
1442      IF( iom_use('snow_ao_cea') )   &
1443         CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:)             )   ! Snow        over ice-free ocean  (cell average)
1444      IF( iom_use('snow_ai_cea') )   &
1445         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
1446
1447      !                                                      ! ========================= !
1448      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1449      !                                                      ! ========================= !
1450      CASE( 'oce only' )                                     ! the required field is directly provided
1451         zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
1452      CASE( 'conservative' )                                      ! the required fields are directly provided
1453         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1454         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1455            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1456         ELSE
1457            ! Set all category values equal for the moment
1458            DO jl=1,jpl
1459               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1460            ENDDO
1461         ENDIF
1462      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
1463         zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1464         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1465            DO jl=1,jpl
1466               zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1467               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1468            ENDDO
1469         ELSE
1470            qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1471            DO jl=1,jpl
1472               zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1473               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1474            ENDDO
1475         ENDIF
1476      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
1477! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1478         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1479         zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1480            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1481            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
1482      END SELECT
1483!!gm
1484!!    currently it is taken into account in leads budget but not in the zqns_tot, and thus not in
1485!!    the flux that enter the ocean....
1486!!    moreover 1 - it is not diagnose anywhere....
1487!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1488!!
1489!! similar job should be done for snow and precipitation temperature
1490      !                                     
1491      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting
1492         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting
1493         zqns_tot(:,:) = zqns_tot(:,:) - ztmp(:,:)
1494         IF( iom_use('hflx_cal_cea') )   &
1495            CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1496      ENDIF
1497
1498      ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus
1499      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1500
1501#if defined key_lim3
1502      CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
1503
1504      ! --- evaporation --- !
1505      ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation
1506      ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice
1507      !                 but it is incoherent WITH the ice model 
1508      DO jl=1,jpl
1509         evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1)
1510      ENDDO
1511      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean
1512
1513      ! --- evaporation minus precipitation --- !
1514      emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:)
1515
1516      ! --- non solar flux over ocean --- !
1517      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
1518      zqns_oce = 0._wp
1519      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
1520
1521      ! --- heat flux associated with emp --- !
1522      zsnw(:,:) = 0._wp
1523      CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing
1524      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
1525         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip
1526         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean
1527      qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
1528         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
1529
1530      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
1531      zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
1532
1533      ! --- total non solar flux --- !
1534      zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:)
1535
1536      ! --- in case both coupled/forced are active, we must mix values --- !
1537      IF( ln_mixcpl ) THEN
1538         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
1539         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
1540         DO jl=1,jpl
1541            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
1542         ENDDO
1543         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
1544         qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
1545!!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0)
1546      ELSE
1547         qns_tot  (:,:  ) = zqns_tot  (:,:  )
1548         qns_oce  (:,:  ) = zqns_oce  (:,:  )
1549         qns_ice  (:,:,:) = zqns_ice  (:,:,:)
1550         qprec_ice(:,:)   = zqprec_ice(:,:)
1551         qemp_oce (:,:)   = zqemp_oce (:,:)
1552      ENDIF
1553
1554      CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
1555#else
1556
1557      ! clem: this formulation is certainly wrong... but better than it was...
1558      zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with:
1559         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
1560         &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
1561         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
1562
1563     IF( ln_mixcpl ) THEN
1564         qns_tot(:,:) = qns(:,:) * p_frld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
1565         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
1566         DO jl=1,jpl
1567            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
1568         ENDDO
1569      ELSE
1570         qns_tot(:,:  ) = zqns_tot(:,:  )
1571         qns_ice(:,:,:) = zqns_ice(:,:,:)
1572      ENDIF
1573
1574#endif
1575
1576      !                                                      ! ========================= !
1577      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1578      !                                                      ! ========================= !
1579      CASE( 'oce only' )
1580         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
1581      CASE( 'conservative' )
1582         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1583         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1584            zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1585         ELSE
1586            ! Set all category values equal for the moment
1587            DO jl=1,jpl
1588               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1589            ENDDO
1590         ENDIF
1591         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1592         zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1593      CASE( 'oce and ice' )
1594         zqsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1595         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1596            DO jl=1,jpl
1597               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1598               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1599            ENDDO
1600         ELSE
1601            qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1602            DO jl=1,jpl
1603               zqsr_tot(:,:   ) = zqsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1604               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1605            ENDDO
1606         ENDIF
1607      CASE( 'mixed oce-ice' )
1608         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1609! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1610!       Create solar heat flux over ice using incoming solar heat flux and albedos
1611!       ( see OASIS3 user guide, 5th edition, p39 )
1612         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1613            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1614            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1615      END SELECT
1616      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
1617         zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
1618         DO jl=1,jpl
1619            zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
1620         ENDDO
1621      ENDIF
1622
1623#if defined key_lim3
1624      CALL wrk_alloc( jpi,jpj, zqsr_oce ) 
1625      ! --- solar flux over ocean --- !
1626      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
1627      zqsr_oce = 0._wp
1628      WHERE( p_frld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / p_frld(:,:)
1629
1630      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
1631      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
1632
1633      CALL wrk_dealloc( jpi,jpj, zqsr_oce ) 
1634#endif
1635
1636      IF( ln_mixcpl ) THEN
1637         qsr_tot(:,:) = qsr(:,:) * p_frld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
1638         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
1639         DO jl=1,jpl
1640            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
1641         ENDDO
1642      ELSE
1643         qsr_tot(:,:  ) = zqsr_tot(:,:  )
1644         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
1645      ENDIF
1646
1647      !                                                      ! ========================= !
1648      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
1649      !                                                      ! ========================= !
1650      CASE ('coupled')
1651         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1652            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
1653         ELSE
1654            ! Set all category values equal for the moment
1655            DO jl=1,jpl
1656               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
1657            ENDDO
1658         ENDIF
1659      END SELECT
1660     
1661      IF( ln_mixcpl ) THEN
1662         DO jl=1,jpl
1663            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
1664         ENDDO
1665      ELSE
1666         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
1667      ENDIF
1668     
1669      !                                                      ! ========================= !
1670      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
1671      !                                                      ! ========================= !
1672      CASE ('coupled')
1673         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
1674         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
1675      END SELECT
1676
1677      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
1678      ! Used for LIM2 and LIM3
1679      ! Coupled case: since cloud cover is not received from atmosphere
1680      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
1681      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
1682      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
1683
1684      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
1685      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
1686      !
1687      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
1688      !
1689   END SUBROUTINE sbc_cpl_ice_flx
1690   
1691   
1692   SUBROUTINE sbc_cpl_snd( kt )
1693      !!----------------------------------------------------------------------
1694      !!             ***  ROUTINE sbc_cpl_snd  ***
1695      !!
1696      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
1697      !!
1698      !! ** Method  :   send to the atmosphere through a call to cpl_snd
1699      !!              all the needed fields (as defined in sbc_cpl_init)
1700      !!----------------------------------------------------------------------
1701      INTEGER, INTENT(in) ::   kt
1702      !
1703      INTEGER ::   ji, jj, jl   ! dummy loop indices
1704      INTEGER ::   isec, info   ! local integer
1705      REAL(wp) ::   zumax, zvmax
1706      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
1707      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
1708      !!----------------------------------------------------------------------
1709      !
1710      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
1711      !
1712      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1713      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1714
1715      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
1716
1717      zfr_l(:,:) = 1.- fr_i(:,:)
1718      !                                                      ! ------------------------- !
1719      !                                                      !    Surface temperature    !   in Kelvin
1720      !                                                      ! ------------------------- !
1721      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
1722         
1723         IF ( nn_components == jp_iam_opa ) THEN
1724            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of ln_useCT on the received part
1725         ELSE
1726            ! we must send the surface potential temperature
1727            IF( ln_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
1728            ELSE                    ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
1729            ENDIF
1730            !
1731            SELECT CASE( sn_snd_temp%cldes)
1732            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
1733            CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
1734               SELECT CASE( sn_snd_temp%clcat )
1735               CASE( 'yes' )   
1736                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl)
1737               CASE( 'no' )
1738                  WHERE( SUM( a_i, dim=3 ) /= 0. )
1739                     ztmp3(:,:,1) = SUM( tn_ice * a_i, dim=3 ) / SUM( a_i, dim=3 )
1740                  ELSEWHERE
1741                     ztmp3(:,:,1) = rt0 ! TODO: Is freezing point a good default? (Maybe SST is better?)
1742                  END WHERE
1743               CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1744               END SELECT
1745            CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)   
1746               SELECT CASE( sn_snd_temp%clcat )
1747               CASE( 'yes' )   
1748                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1749               CASE( 'no' )
1750                  ztmp3(:,:,:) = 0.0
1751                  DO jl=1,jpl
1752                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
1753                  ENDDO
1754               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1755               END SELECT
1756            CASE( 'mixed oce-ice'        )   
1757               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) 
1758               DO jl=1,jpl
1759                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
1760               ENDDO
1761            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
1762            END SELECT
1763         ENDIF
1764         IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1765         IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
1766         IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1767      ENDIF
1768      !                                                      ! ------------------------- !
1769      !                                                      !           Albedo          !
1770      !                                                      ! ------------------------- !
1771      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
1772         SELECT CASE( sn_snd_alb%cldes )
1773         CASE( 'ice'          )   ; ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
1774         CASE( 'weighted ice' )   ; ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1775         CASE default             ; CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
1776         END SELECT
1777         CALL cpl_snd( jps_albice, isec, ztmp3, info )
1778      ENDIF
1779      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
1780         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
1781         DO jl=1,jpl
1782            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
1783         ENDDO
1784         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1785      ENDIF
1786      !                                                      ! ------------------------- !
1787      !                                                      !  Ice fraction & Thickness !
1788      !                                                      ! ------------------------- !
1789      ! Send ice fraction field to atmosphere
1790      IF( ssnd(jps_fice)%laction ) THEN
1791         SELECT CASE( sn_snd_thick%clcat )
1792         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
1793         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
1794         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1795         END SELECT
1796         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
1797      ENDIF
1798     
1799      ! Send ice fraction field to OPA (sent by SAS in SAS-OPA coupling)
1800      IF( ssnd(jps_fice2)%laction ) THEN
1801         ztmp3(:,:,1) = fr_i(:,:)
1802         IF( ssnd(jps_fice2)%laction )   CALL cpl_snd( jps_fice2, isec, ztmp3, info )
1803      ENDIF
1804
1805      ! Send ice and snow thickness field
1806      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
1807         SELECT CASE( sn_snd_thick%cldes)
1808         CASE( 'none'                  )       ! nothing to do
1809         CASE( 'weighted ice and snow' )   
1810            SELECT CASE( sn_snd_thick%clcat )
1811            CASE( 'yes' )   
1812               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
1813               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
1814            CASE( 'no' )
1815               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
1816               DO jl=1,jpl
1817                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
1818                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
1819               ENDDO
1820            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1821            END SELECT
1822         CASE( 'ice and snow'         )   
1823            SELECT CASE( sn_snd_thick%clcat )
1824            CASE( 'yes' )
1825               ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
1826               ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
1827            CASE( 'no' )
1828               WHERE( SUM( a_i, dim=3 ) /= 0. )
1829                  ztmp3(:,:,1) = SUM( ht_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
1830                  ztmp4(:,:,1) = SUM( ht_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
1831               ELSEWHERE
1832                 ztmp3(:,:,1) = 0.
1833                 ztmp4(:,:,1) = 0.
1834               END WHERE
1835            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1836            END SELECT
1837         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
1838         END SELECT
1839         IF( ssnd(jps_hice)%laction )   CALL cpl_snd( jps_hice, isec, ztmp3, info )
1840         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
1841      ENDIF
1842      !
1843#if defined key_cpl_carbon_cycle
1844      !                                                      ! ------------------------- !
1845      !                                                      !  CO2 flux from PISCES     !
1846      !                                                      ! ------------------------- !
1847      IF( ssnd(jps_co2)%laction )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
1848      !
1849#endif
1850      !                                                      ! ------------------------- !
1851      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
1852         !                                                   ! ------------------------- !
1853         !   
1854         !                                                  j+1   j     -----V---F
1855         ! surface velocity always sent from T point                     !       |
1856         !                                                        j      |   T   U
1857         !                                                               |       |
1858         !                                                   j    j-1   -I-------|
1859         !                                               (for I)         |       |
1860         !                                                              i-1  i   i
1861         !                                                               i      i+1 (for I)
1862         IF( nn_components == jp_iam_opa ) THEN
1863            zotx1(:,:) = un(:,:,1) 
1864            zoty1(:,:) = vn(:,:,1) 
1865         ELSE       
1866            SELECT CASE( TRIM( sn_snd_crt%cldes ) )
1867            CASE( 'oce only'             )      ! C-grid ==> T
1868               DO jj = 2, jpjm1
1869                  DO ji = fs_2, fs_jpim1   ! vector opt.
1870                     zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
1871                     zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
1872                  END DO
1873               END DO
1874            CASE( 'weighted oce and ice' )   
1875               SELECT CASE ( cp_ice_msh )
1876               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1877                  DO jj = 2, jpjm1
1878                     DO ji = fs_2, fs_jpim1   ! vector opt.
1879                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1880                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
1881                        zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1882                        zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1883                     END DO
1884                  END DO
1885               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1886                  DO jj = 2, jpjm1
1887                     DO ji = 2, jpim1   ! NO vector opt.
1888                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1889                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1890                        zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1891                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1892                        zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1893                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1894                     END DO
1895                  END DO
1896               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1897                  DO jj = 2, jpjm1
1898                     DO ji = 2, jpim1   ! NO vector opt.
1899                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1900                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1901                        zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1902                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1903                        zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1904                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1905                     END DO
1906                  END DO
1907               END SELECT
1908               CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
1909            CASE( 'mixed oce-ice'        )
1910               SELECT CASE ( cp_ice_msh )
1911               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1912                  DO jj = 2, jpjm1
1913                     DO ji = fs_2, fs_jpim1   ! vector opt.
1914                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
1915                           &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1916                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
1917                           &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1918                     END DO
1919                  END DO
1920               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1921                  DO jj = 2, jpjm1
1922                     DO ji = 2, jpim1   ! NO vector opt.
1923                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1924                           &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1925                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1926                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1927                           &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1928                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1929                     END DO
1930                  END DO
1931               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1932                  DO jj = 2, jpjm1
1933                     DO ji = 2, jpim1   ! NO vector opt.
1934                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1935                           &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1936                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1937                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1938                           &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1939                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1940                     END DO
1941                  END DO
1942               END SELECT
1943            END SELECT
1944            CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
1945            !
1946         ENDIF
1947         !
1948         !
1949         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
1950            !                                                                     ! Ocean component
1951            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
1952            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
1953            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
1954            zoty1(:,:) = ztmp2(:,:)
1955            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
1956               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
1957               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
1958               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
1959               zity1(:,:) = ztmp2(:,:)
1960            ENDIF
1961         ENDIF
1962         !
1963         ! spherical coordinates to cartesian -> 2 components to 3 components
1964         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
1965            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
1966            ztmp2(:,:) = zoty1(:,:)
1967            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
1968            !
1969            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
1970               ztmp1(:,:) = zitx1(:,:)
1971               ztmp1(:,:) = zity1(:,:)
1972               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
1973            ENDIF
1974         ENDIF
1975         !
1976         IF( ssnd(jps_ocx1)%laction )   CALL cpl_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
1977         IF( ssnd(jps_ocy1)%laction )   CALL cpl_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
1978         IF( ssnd(jps_ocz1)%laction )   CALL cpl_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
1979         !
1980         IF( ssnd(jps_ivx1)%laction )   CALL cpl_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
1981         IF( ssnd(jps_ivy1)%laction )   CALL cpl_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
1982         IF( ssnd(jps_ivz1)%laction )   CALL cpl_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
1983         !
1984      ENDIF
1985      !
1986      !
1987      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
1988      !                                                        ! SSH
1989      IF( ssnd(jps_ssh )%laction )  THEN
1990         !                          ! removed inverse barometer ssh when Patm
1991         !                          forcing is used (for sea-ice dynamics)
1992         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
1993         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
1994         ENDIF
1995         CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
1996
1997      ENDIF
1998      !                                                        ! SSS
1999      IF( ssnd(jps_soce  )%laction )  THEN
2000         CALL cpl_snd( jps_soce  , isec, RESHAPE ( tsn(:,:,1,jp_sal), (/jpi,jpj,1/) ), info )
2001      ENDIF
2002      !                                                        ! first T level thickness
2003      IF( ssnd(jps_e3t1st )%laction )  THEN
2004         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( fse3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
2005      ENDIF
2006      !                                                        ! Qsr fraction
2007      IF( ssnd(jps_fraqsr)%laction )  THEN
2008         CALL cpl_snd( jps_fraqsr, isec, RESHAPE ( fraqsr_1lev(:,:) , (/jpi,jpj,1/) ), info )
2009      ENDIF
2010      !
2011      !  Fields sent by SAS to OPA when OASIS coupling
2012      !                                                        ! Solar heat flux
2013      IF( ssnd(jps_qsroce)%laction )  CALL cpl_snd( jps_qsroce, isec, RESHAPE ( qsr , (/jpi,jpj,1/) ), info )
2014      IF( ssnd(jps_qnsoce)%laction )  CALL cpl_snd( jps_qnsoce, isec, RESHAPE ( qns , (/jpi,jpj,1/) ), info )
2015      IF( ssnd(jps_oemp  )%laction )  CALL cpl_snd( jps_oemp  , isec, RESHAPE ( emp , (/jpi,jpj,1/) ), info )
2016      IF( ssnd(jps_sflx  )%laction )  CALL cpl_snd( jps_sflx  , isec, RESHAPE ( sfx , (/jpi,jpj,1/) ), info )
2017      IF( ssnd(jps_otx1  )%laction )  CALL cpl_snd( jps_otx1  , isec, RESHAPE ( utau, (/jpi,jpj,1/) ), info )
2018      IF( ssnd(jps_oty1  )%laction )  CALL cpl_snd( jps_oty1  , isec, RESHAPE ( vtau, (/jpi,jpj,1/) ), info )
2019      IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info )
2020      IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
2021
2022      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
2023      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
2024      !
2025      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
2026      !
2027   END SUBROUTINE sbc_cpl_snd
2028   
2029   !!======================================================================
2030END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.