source: branches/UKMO/dev_r5518_rm_um_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5856

Last change on this file since 5856 was 5856, checked in by jcastill, 5 years ago

The Surface stress wind is only received from atmosphere if specified in namsbc_cpl

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