New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbccpl.F90 in branches/UKMO/dev_r5518_CICE_coupling_GSI7/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_CICE_coupling_GSI7/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5612

Last change on this file since 5612 was 5612, checked in by dancopsey, 9 years ago

Merged with revision 5518 of the trunk (NEMO3.6_stable).

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