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/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 4733

Last change on this file since 4733 was 4733, checked in by vancop, 10 years ago

Fix energy budget in coupled case

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