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

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

#1260 : see the associated wiki page for explanation

  • Property svn:keywords set to Id
File size: 95.7 KB
Line 
1MODULE sbccpl
2   !!======================================================================
3   !!                       ***  MODULE  sbccpl  ***
4   !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode
5   !!======================================================================
6   !! History :  2.0  ! 2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
7   !!            3.0  ! 2008-02  (G. Madec, C Talandier)  surface module
8   !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
9   !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
10   !!----------------------------------------------------------------------
11#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 recieved information from
213      !!                the atmospheric component
214      !!
215      !! ** Method  : * Read namsbc_cpl namelist
216      !!              * define the receive interface
217      !!              * define the send    interface
218      !!              * initialise the OASIS coupler
219      !!----------------------------------------------------------------------
220      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
221      !!
222      INTEGER ::   jn   ! dummy loop index
223      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      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            DO jj = 2, jpjm1
755               DO ji = fs_2, fs_jpim1   ! vect. opt.
756                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
757                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
758                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
759               END DO
760            END DO
761            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
762            llnewtau = .TRUE.
763         ELSE
764            llnewtau = .FALSE.
765         ENDIF
766      ELSE
767         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
768         ! Stress module can be negative when received (interpolation problem)
769         IF( llnewtau ) THEN
770            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
771         ENDIF
772      ENDIF
773     
774      !                                                      ! ========================= !
775      !                                                      !      10 m wind speed      !   (wndm)
776      !                                                      ! ========================= !
777      !
778      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
779         ! => need to be done only when taumod was changed
780         IF( llnewtau ) THEN
781            zcoef = 1. / ( zrhoa * zcdrag ) 
782            DO jj = 1, jpj
783               DO ji = 1, jpi 
784                  wndm(ji,jj) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
785               END DO
786            END DO
787         ENDIF
788      ELSE
789         IF ( nrcvinfo(jpr_w10m) == OASIS_Rcv ) wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
790      ENDIF
791
792      ! u(v)tau and taum will be modified by ice model
793      ! -> need to be reset before each call of the ice/fsbc     
794      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
795         utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
796         vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
797         taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
798         CALL iom_put( "taum_oce", taum )   ! output wind stress module
799      ENDIF
800
801#if defined key_cpl_carbon_cycle
802      !                                                              ! atmosph. CO2 (ppm)
803      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
804#endif
805
806      !                                                      ! ========================= !
807      IF( k_ice <= 1 ) THEN                                  !  heat & freshwater fluxes ! (Ocean only case)
808         !                                                   ! ========================= !
809         !
810         !                                                       ! total freshwater fluxes over the ocean (emp)
811         SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
812         CASE( 'conservative' )
813            emp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
814         CASE( 'oce only', 'oce and ice' )
815            emp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
816         CASE default
817            CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
818         END SELECT
819         !
820         !                                                        ! runoffs and calving (added in emp)
821         IF( srcv(jpr_rnf)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
822         IF( srcv(jpr_cal)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_cal)%z3(:,:,1)
823         !
824!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
825!!gm                                       at least should be optional...
826!!         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN     ! add to the total freshwater budget
827!!            ! remove negative runoff
828!!            zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1e2t(:,:) * tmask_i(:,:) )
829!!            zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1e2t(:,:) * tmask_i(:,:) )
830!!            IF( lk_mpp )   CALL mpp_sum( zcumulpos )   ! sum over the global domain
831!!            IF( lk_mpp )   CALL mpp_sum( zcumulneg )
832!!            IF( zcumulpos /= 0. ) THEN                 ! distribute negative runoff on positive runoff grid points
833!!               zcumulneg = 1.e0 + zcumulneg / zcumulpos
834!!               frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
835!!            ENDIF     
836!!            ! add runoff to e-p
837!!            emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
838!!         ENDIF
839!!gm  end of internal cooking
840         !
841         !                                                       ! non solar heat flux over the ocean (qns)
842         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
843         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
844         ! add the latent heat of solid precip. melting
845         IF( srcv(jpr_snow  )%laction )   THEN                         ! update qns over the free ocean with:
846              qns(:,:) = qns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus  & ! energy for melting solid precipitation over the free ocean
847           &           - emp(:,:) * sst_m(:,:) * rcp                   ! remove heat content due to mass flux (assumed to be at SST)
848         ENDIF
849
850         !                                                       ! solar flux over the ocean          (qsr)
851         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
852         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
853         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
854         !
855 
856      ENDIF
857      !
858      CALL wrk_dealloc( jpi,jpj, ztx, zty )
859      !
860      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
861      !
862   END SUBROUTINE sbc_cpl_rcv
863   
864
865   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
866      !!----------------------------------------------------------------------
867      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
868      !!
869      !! ** Purpose :   provide the stress over sea-ice in coupled mode
870      !!
871      !! ** Method  :   transform the received stress from the atmosphere into
872      !!             an atmosphere-ice stress in the (i,j) ocean referencial
873      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
874      !!                'C'-grid : i- (j-) components given at U- (V-) point
875      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
876      !!
877      !!                The received stress are :
878      !!                 - defined by 3 components (if cartesian coordinate)
879      !!                        or by 2 components (if spherical)
880      !!                 - oriented along geographical   coordinate (if eastward-northward)
881      !!                        or  along the local grid coordinate (if local grid)
882      !!                 - given at U- and V-point, resp.   if received on 2 grids
883      !!                        or at a same point (T or I) if received on 1 grid
884      !!                Therefore and if necessary, they are successively
885      !!             processed in order to obtain them
886      !!                 first  as  2 components on the sphere
887      !!                 second as  2 components oriented along the local grid
888      !!                 third  as  2 components on the cp_ice_msh point
889      !!
890      !!                Except in 'oce and ice' case, only one vector stress field
891      !!             is received. It has already been processed in sbc_cpl_rcv
892      !!             so that it is now defined as (i,j) components given at U-
893      !!             and V-points, respectively. Therefore, only the third
894      !!             transformation is done and only if the ice-grid is a 'I'-grid.
895      !!
896      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
897      !!----------------------------------------------------------------------
898      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
899      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
900      !!
901      INTEGER ::   ji, jj                          ! dummy loop indices
902      INTEGER ::   itx                             ! index of taux over ice
903      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
904      !!----------------------------------------------------------------------
905      !
906      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
907      !
908      CALL wrk_alloc( jpi,jpj, ztx, zty )
909
910!AC Pour eviter un stress nul sur la glace dans le cas mixed oce-ice
911      IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN   ;   itx =  jpr_itx1   
912      ELSE                                ;   itx =  jpr_otx1
913      ENDIF
914
915      ! do something only if we just received the stress from atmosphere
916      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
917
918         !                                                                                              ! ======================= !
919!AC Pour eviter un stress nul sur la glace dans le cas mixes oce-ice
920         IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN               !   ice stress received   !
921            !                                                                                           ! ======================= !
922           
923            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
924               !                                                       ! (cartesian to spherical -> 3 to 2 components)
925               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
926                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
927               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
928               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
929               !
930               IF( srcv(jpr_itx2)%laction ) THEN
931                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
932                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
933                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
934                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
935               ENDIF
936               !
937            ENDIF
938            !
939            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
940               !                                                       ! (geographical to local grid -> rotate the components)
941               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
942               IF( srcv(jpr_itx2)%laction ) THEN
943                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
944               ELSE
945                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
946               ENDIF
947               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
948               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
949            ENDIF
950            !                                                   ! ======================= !
951         ELSE                                                   !     use ocean stress    !
952            !                                                   ! ======================= !
953            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
954            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
955            !
956         ENDIF
957
958         !                                                      ! ======================= !
959         !                                                      !     put on ice grid     !
960         !                                                      ! ======================= !
961         !   
962         !                                                  j+1   j     -----V---F
963         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
964         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
965         !                                                               |       |
966         !                                                   j    j-1   -I-------|
967         !                                               (for I)         |       |
968         !                                                              i-1  i   i
969         !                                                               i      i+1 (for I)
970         SELECT CASE ( cp_ice_msh )
971            !
972         CASE( 'I' )                                         ! B-grid ==> I
973            SELECT CASE ( srcv(jpr_itx1)%clgrid )
974            CASE( 'U' )
975               DO jj = 2, jpjm1                                   ! (U,V) ==> I
976                  DO ji = 2, jpim1   ! NO vector opt.
977                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
978                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
979                  END DO
980               END DO
981            CASE( 'F' )
982               DO jj = 2, jpjm1                                   ! F ==> I
983                  DO ji = 2, jpim1   ! NO vector opt.
984                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
985                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
986                  END DO
987               END DO
988            CASE( 'T' )
989               DO jj = 2, jpjm1                                   ! T ==> I
990                  DO ji = 2, jpim1   ! NO vector opt.
991                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
992                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
993                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
994                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
995                  END DO
996               END DO
997            CASE( 'I' )
998               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
999               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1000            END SELECT
1001            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1002               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1003            ENDIF
1004            !
1005         CASE( 'F' )                                         ! B-grid ==> F
1006            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1007            CASE( 'U' )
1008               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1009                  DO ji = fs_2, fs_jpim1   ! vector opt.
1010                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1011                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1012                  END DO
1013               END DO
1014            CASE( 'I' )
1015               DO jj = 2, jpjm1                                   ! I ==> F
1016                  DO ji = 2, jpim1   ! NO vector opt.
1017                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1018                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1019                  END DO
1020               END DO
1021            CASE( 'T' )
1022               DO jj = 2, jpjm1                                   ! T ==> F
1023                  DO ji = 2, jpim1   ! NO vector opt.
1024                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1025                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1026                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1027                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1028                  END DO
1029               END DO
1030            CASE( 'F' )
1031               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1032               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1033            END SELECT
1034            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1035               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1036            ENDIF
1037            !
1038         CASE( 'C' )                                         ! C-grid ==> U,V
1039            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1040            CASE( 'U' )
1041               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1042               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1043            CASE( 'F' )
1044               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1045                  DO ji = fs_2, fs_jpim1   ! vector opt.
1046                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1047                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1048                  END DO
1049               END DO
1050            CASE( 'T' )
1051               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1052                  DO ji = fs_2, fs_jpim1   ! vector opt.
1053                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1054                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1055                  END DO
1056               END DO
1057            CASE( 'I' )
1058               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1059                  DO ji = 2, jpim1   ! NO vector opt.
1060                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1061                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1062                  END DO
1063               END DO
1064            END SELECT
1065            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1066               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1067            ENDIF
1068         END SELECT
1069
1070      ENDIF
1071      !   
1072      CALL wrk_dealloc( jpi,jpj, ztx, zty )
1073      !
1074      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
1075      !
1076   END SUBROUTINE sbc_cpl_ice_tau
1077   
1078
1079   SUBROUTINE sbc_cpl_ice_flx( p_frld  , palbi   , psst    , pist    )
1080      !!----------------------------------------------------------------------
1081      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1082      !!
1083      !! ** Purpose :   provide the heat and freshwater fluxes of the
1084      !!              ocean-ice system.
1085      !!
1086      !! ** Method  :   transform the fields received from the atmosphere into
1087      !!             surface heat and fresh water boundary condition for the
1088      !!             ice-ocean system. The following fields are provided:
1089      !!              * total non solar, solar and freshwater fluxes (qns_tot,
1090      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1091      !!             NB: emp_tot include runoffs and calving.
1092      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1093      !!             emp_ice = sublimation - solid precipitation as liquid
1094      !!             precipitation are re-routed directly to the ocean and
1095      !!             runoffs and calving directly enter the ocean.
1096      !!              * solid precipitation (sprecip), used to add to qns_tot
1097      !!             the heat lost associated to melting solid precipitation
1098      !!             over the ocean fraction.
1099      !!       ===>> CAUTION here this changes the net heat flux received from
1100      !!             the atmosphere
1101      !!
1102      !!                  - the fluxes have been separated from the stress as
1103      !!                 (a) they are updated at each ice time step compare to
1104      !!                 an update at each coupled time step for the stress, and
1105      !!                 (b) the conservative computation of the fluxes over the
1106      !!                 sea-ice area requires the knowledge of the ice fraction
1107      !!                 after the ice advection and before the ice thermodynamics,
1108      !!                 so that the stress is updated before the ice dynamics
1109      !!                 while the fluxes are updated after it.
1110      !!
1111      !! ** Action  :   update at each nf_ice time step:
1112      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1113      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1114      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
1115      !!                   emp_ice            ice sublimation - solid precipitation over the ice
1116      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
1117      !!                   sprecip             solid precipitation over the ocean 
1118      !!----------------------------------------------------------------------
1119      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1120      ! optional arguments, used only in 'mixed oce-ice' case
1121      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1122      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celcius]
1123      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature     [Kelvin]
1124      !
1125      INTEGER ::   jl   ! dummy loop index
1126      REAL(wp), POINTER, DIMENSION(:,:) ::   zcptn, ztmp, zicefr
1127      !!----------------------------------------------------------------------
1128      !
1129      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1130      !
1131      CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr )
1132
1133      zicefr(:,:) = 1.- p_frld(:,:)
1134      zcptn(:,:) = rcp * sst_m(:,:)
1135      !
1136      !                                                      ! ========================= !
1137      !                                                      !    freshwater budget      !   (emp)
1138      !                                                      ! ========================= !
1139      !
1140      !                                                           ! total Precipitations - total Evaporation (emp_tot)
1141      !                                                           ! solid precipitation  - sublimation       (emp_ice)
1142      !                                                           ! solid Precipitation                      (sprecip)
1143      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1144      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1145         sprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                 ! May need to ensure positive here
1146         tprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + sprecip (:,:) ! May need to ensure positive here
1147         emp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - tprecip(:,:)
1148         emp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1149                           CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation
1150         IF( lk_diaar5 )   CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
1151         ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1152                           CALL iom_put( 'evap_ao_cea'  , ztmp                            )   ! ice-free oce evap (cell average)
1153         IF( lk_diaar5 )   CALL iom_put( 'hflx_evap_cea', ztmp(:,:         ) * zcptn(:,:) )   ! heat flux from from evap (cell ave)
1154      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1155         emp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1156         emp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
1157         sprecip(:,:) = - frcv(jpr_semp)%z3(:,:,1) + frcv(jpr_ievp)%z3(:,:,1)
1158      END SELECT
1159
1160      CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
1161      CALL iom_put( 'snow_ao_cea', sprecip(:,:         ) * p_frld(:,:)    )   ! Snow        over ice-free ocean  (cell average)
1162      CALL iom_put( 'snow_ai_cea', sprecip(:,:         ) * zicefr(:,:)    )   ! Snow        over sea-ice         (cell average)
1163      CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
1164      !   
1165      !                                                           ! runoffs and calving (put in emp_tot)
1166      IF( srcv(jpr_rnf)%laction ) THEN
1167         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)
1168                           CALL iom_put( 'runoffs'      , frcv(jpr_rnf)%z3(:,:,1)              )   ! rivers
1169         IF( lk_diaar5 )   CALL iom_put( 'hflx_rnf_cea' , frcv(jpr_rnf)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from rivers
1170      ENDIF
1171      IF( srcv(jpr_cal)%laction ) THEN
1172         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1173         CALL iom_put( 'calving', frcv(jpr_cal)%z3(:,:,1) )
1174      ENDIF
1175      !
1176!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
1177!!gm                                       at least should be optional...
1178!!       ! remove negative runoff                            ! sum over the global domain
1179!!       zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1e2t(:,:) * tmask_i(:,:) )
1180!!       zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1e2t(:,:) * tmask_i(:,:) )
1181!!       IF( lk_mpp )   CALL mpp_sum( zcumulpos )
1182!!       IF( lk_mpp )   CALL mpp_sum( zcumulneg )
1183!!       IF( zcumulpos /= 0. ) THEN                          ! distribute negative runoff on positive runoff grid points
1184!!          zcumulneg = 1.e0 + zcumulneg / zcumulpos
1185!!          frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
1186!!       ENDIF     
1187!!       emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)   ! add runoff to e-p
1188!!
1189!!gm  end of internal cooking
1190
1191      !                                                      ! ========================= !
1192      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1193      !                                                      ! ========================= !
1194      CASE( 'oce only' )                                     ! the required field is directly provided
1195         qns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
1196      CASE( 'conservative' )                                      ! the required fields are directly provided
1197         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1198         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1199            qns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1200         ELSE
1201            ! Set all category values equal for the moment
1202            DO jl=1,jpl
1203               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1204            ENDDO
1205         ENDIF
1206      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
1207         qns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1208         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1209            DO jl=1,jpl
1210               qns_tot(:,:   ) = qns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1211               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1212            ENDDO
1213         ELSE
1214            DO jl=1,jpl
1215               qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1216               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1217            ENDDO
1218         ENDIF
1219      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
1220! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1221         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1222         qns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1223            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1224            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
1225      END SELECT
1226      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus
1227      qns_tot(:,:) = qns_tot(:,:)                         &            ! qns_tot update over free ocean with:
1228         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
1229         &          - (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
1230         &             - emp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
1231      IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1232!!gm
1233!!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in
1234!!    the flux that enter the ocean....
1235!!    moreover 1 - it is not diagnose anywhere....
1236!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1237!!
1238!! similar job should be done for snow and precipitation temperature
1239      !                                     
1240      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting
1241         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting
1242         qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)
1243         IF( lk_diaar5 )   CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1244      ENDIF
1245
1246      !                                                      ! ========================= !
1247      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1248      !                                                      ! ========================= !
1249      CASE( 'oce only' )
1250         qsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
1251      CASE( 'conservative' )
1252         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1253         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1254            qsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1255         ELSE
1256            ! Set all category values equal for the moment
1257            DO jl=1,jpl
1258               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1259            ENDDO
1260         ENDIF
1261         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1262         qsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1263      CASE( 'oce and ice' )
1264         qsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1265         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1266            DO jl=1,jpl
1267               qsr_tot(:,:   ) = qsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1268               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1269            ENDDO
1270         ELSE
1271            DO jl=1,jpl
1272               qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1273               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1274            ENDDO
1275         ENDIF
1276      CASE( 'mixed oce-ice' )
1277         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1278! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1279!       Create solar heat flux over ice using incoming solar heat flux and albedos
1280!       ( see OASIS3 user guide, 5th edition, p39 )
1281         qsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1282            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1283            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1284      END SELECT
1285      IF( ln_dm2dc ) THEN   ! modify qsr to include the diurnal cycle
1286         qsr_tot(:,:  ) = sbc_dcy( qsr_tot(:,:  ) )
1287         DO jl=1,jpl
1288            qsr_ice(:,:,jl) = sbc_dcy( qsr_ice(:,:,jl) )
1289         ENDDO
1290      ENDIF
1291
1292      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )
1293      CASE ('coupled')
1294         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1295            dqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
1296         ELSE
1297            ! Set all category values equal for the moment
1298            DO jl=1,jpl
1299               dqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
1300            ENDDO
1301         ENDIF
1302      END SELECT
1303
1304      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )
1305      CASE ('coupled')
1306         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
1307         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
1308      END SELECT
1309
1310      !    Ice Qsr penetration used (only?)in lim2 or lim3
1311      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer
1312      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
1313      ! Coupled case: since cloud cover is not received from atmosphere
1314      !               ===> defined as constant value -> definition done in sbc_cpl_init
1315      fr1_i0(:,:) = 0.18
1316      fr2_i0(:,:) = 0.82
1317
1318
1319      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
1320      !
1321      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
1322      !
1323   END SUBROUTINE sbc_cpl_ice_flx
1324   
1325   
1326   SUBROUTINE sbc_cpl_snd( kt )
1327      !!----------------------------------------------------------------------
1328      !!             ***  ROUTINE sbc_cpl_snd  ***
1329      !!
1330      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
1331      !!
1332      !! ** Method  :   send to the atmosphere through a call to cpl_prism_snd
1333      !!              all the needed fields (as defined in sbc_cpl_init)
1334      !!----------------------------------------------------------------------
1335      INTEGER, INTENT(in) ::   kt
1336      !
1337      INTEGER ::   ji, jj, jl   ! dummy loop indices
1338      INTEGER ::   isec, info   ! local integer
1339      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
1340      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
1341      !!----------------------------------------------------------------------
1342      !
1343      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
1344      !
1345      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1346      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1347
1348      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
1349
1350      zfr_l(:,:) = 1.- fr_i(:,:)
1351
1352      !                                                      ! ------------------------- !
1353      !                                                      !    Surface temperature    !   in Kelvin
1354      !                                                      ! ------------------------- !
1355      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
1356         SELECT CASE( sn_snd_temp%cldes)
1357         CASE( 'oce only'             )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0
1358         CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:)   
1359            SELECT CASE( sn_snd_temp%clcat )
1360            CASE( 'yes' )   
1361               ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1362            CASE( 'no' )
1363               ztmp3(:,:,:) = 0.0
1364               DO jl=1,jpl
1365                  ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
1366               ENDDO
1367            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1368            END SELECT
1369         CASE( 'mixed oce-ice'        )   
1370            ztmp1(:,:) = ( tsn(:,:,1,1) + rt0 ) * zfr_l(:,:) 
1371            DO jl=1,jpl
1372               ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
1373            ENDDO
1374         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
1375         END SELECT
1376         IF( ssnd(jps_toce)%laction )   CALL cpl_prism_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1377         IF( ssnd(jps_tice)%laction )   CALL cpl_prism_snd( jps_tice, isec, ztmp3, info )
1378         IF( ssnd(jps_tmix)%laction )   CALL cpl_prism_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1379      ENDIF
1380      !
1381      !                                                      ! ------------------------- !
1382      !                                                      !           Albedo          !
1383      !                                                      ! ------------------------- !
1384      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
1385         ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1386         CALL cpl_prism_snd( jps_albice, isec, ztmp3, info )
1387      ENDIF
1388      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
1389         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
1390         DO jl=1,jpl
1391            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
1392         ENDDO
1393         CALL cpl_prism_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1394      ENDIF
1395      !                                                      ! ------------------------- !
1396      !                                                      !  Ice fraction & Thickness !
1397      !                                                      ! ------------------------- !
1398      ! Send ice fraction field
1399      IF( ssnd(jps_fice)%laction ) THEN
1400         SELECT CASE( sn_snd_thick%clcat )
1401         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
1402         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
1403         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1404         END SELECT
1405         CALL cpl_prism_snd( jps_fice, isec, ztmp3, info )
1406      ENDIF
1407
1408      ! Send ice and snow thickness field
1409      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
1410         SELECT CASE( sn_snd_thick%cldes)
1411         CASE( 'none'                  )       ! nothing to do
1412         CASE( 'weighted ice and snow' )   
1413            SELECT CASE( sn_snd_thick%clcat )
1414            CASE( 'yes' )   
1415               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
1416               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
1417            CASE( 'no' )
1418               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
1419               DO jl=1,jpl
1420                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
1421                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
1422               ENDDO
1423            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1424            END SELECT
1425         CASE( 'ice and snow'         )   
1426            ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
1427            ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
1428         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
1429         END SELECT
1430         IF( ssnd(jps_hice)%laction )   CALL cpl_prism_snd( jps_hice, isec, ztmp3, info )
1431         IF( ssnd(jps_hsnw)%laction )   CALL cpl_prism_snd( jps_hsnw, isec, ztmp4, info )
1432      ENDIF
1433      !
1434#if defined key_cpl_carbon_cycle
1435      !                                                      ! ------------------------- !
1436      !                                                      !  CO2 flux from PISCES     !
1437      !                                                      ! ------------------------- !
1438      IF( ssnd(jps_co2)%laction )   CALL cpl_prism_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
1439      !
1440#endif
1441      !                                                      ! ------------------------- !
1442      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
1443         !                                                   ! ------------------------- !
1444         !   
1445         !                                                  j+1   j     -----V---F
1446         ! surface velocity always sent from T point                     !       |
1447         !                                                        j      |   T   U
1448         !                                                               |       |
1449         !                                                   j    j-1   -I-------|
1450         !                                               (for I)         |       |
1451         !                                                              i-1  i   i
1452         !                                                               i      i+1 (for I)
1453         SELECT CASE( TRIM( sn_snd_crt%cldes ) )
1454         CASE( 'oce only'             )      ! C-grid ==> T
1455            DO jj = 2, jpjm1
1456               DO ji = fs_2, fs_jpim1   ! vector opt.
1457                  zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
1458                  zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
1459               END DO
1460            END DO
1461         CASE( 'weighted oce and ice' )   
1462            SELECT CASE ( cp_ice_msh )
1463            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1464               DO jj = 2, jpjm1
1465                  DO ji = fs_2, fs_jpim1   ! vector opt.
1466                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1467                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
1468                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1469                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1470                  END DO
1471               END DO
1472            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1473               DO jj = 2, jpjm1
1474                  DO ji = 2, jpim1   ! NO vector opt.
1475                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1476                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1477                     zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1478                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1479                     zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1480                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1481                  END DO
1482               END DO
1483            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1484               DO jj = 2, jpjm1
1485                  DO ji = 2, jpim1   ! NO vector opt.
1486                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1487                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1488                     zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1489                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1490                     zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1491                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1492                  END DO
1493               END DO
1494            END SELECT
1495            CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
1496         CASE( 'mixed oce-ice'        )
1497            SELECT CASE ( cp_ice_msh )
1498            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1499               DO jj = 2, jpjm1
1500                  DO ji = fs_2, fs_jpim1   ! vector opt.
1501                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
1502                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1503                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
1504                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1505                  END DO
1506               END DO
1507            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1508               DO jj = 2, jpjm1
1509                  DO ji = 2, jpim1   ! NO vector opt.
1510                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1511                        &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1512                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1513                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1514                        &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1515                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1516                  END DO
1517               END DO
1518            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1519               DO jj = 2, jpjm1
1520                  DO ji = 2, jpim1   ! NO vector opt.
1521                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1522                        &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1523                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1524                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1525                        &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1526                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1527                  END DO
1528               END DO
1529            END SELECT
1530         END SELECT
1531         CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
1532         !
1533         !
1534         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
1535            !                                                                     ! Ocean component
1536            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
1537            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
1538            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
1539            zoty1(:,:) = ztmp2(:,:)
1540            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
1541               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
1542               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
1543               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
1544               zity1(:,:) = ztmp2(:,:)
1545            ENDIF
1546         ENDIF
1547         !
1548         ! spherical coordinates to cartesian -> 2 components to 3 components
1549         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
1550            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
1551            ztmp2(:,:) = zoty1(:,:)
1552            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
1553            !
1554            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
1555               ztmp1(:,:) = zitx1(:,:)
1556               ztmp1(:,:) = zity1(:,:)
1557               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
1558            ENDIF
1559         ENDIF
1560         !
1561         IF( ssnd(jps_ocx1)%laction )   CALL cpl_prism_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
1562         IF( ssnd(jps_ocy1)%laction )   CALL cpl_prism_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
1563         IF( ssnd(jps_ocz1)%laction )   CALL cpl_prism_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
1564         !
1565         IF( ssnd(jps_ivx1)%laction )   CALL cpl_prism_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
1566         IF( ssnd(jps_ivy1)%laction )   CALL cpl_prism_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
1567         IF( ssnd(jps_ivz1)%laction )   CALL cpl_prism_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
1568         !
1569      ENDIF
1570      !
1571      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1572      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1573      !
1574      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
1575      !
1576   END SUBROUTINE sbc_cpl_snd
1577   
1578#else
1579   !!----------------------------------------------------------------------
1580   !!   Dummy module                                            NO coupling
1581   !!----------------------------------------------------------------------
1582   USE par_kind        ! kind definition
1583CONTAINS
1584   SUBROUTINE sbc_cpl_snd( kt )
1585      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt
1586   END SUBROUTINE sbc_cpl_snd
1587   !
1588   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
1589      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt, k_fsbc, k_ice
1590   END SUBROUTINE sbc_cpl_rcv
1591   !
1592   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1593      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1594      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1595      p_taui(:,:) = 0.   ;   p_tauj(:,:) = 0. ! stupid definition to avoid warning message when compiling...
1596      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?'
1597   END SUBROUTINE sbc_cpl_ice_tau
1598   !
1599   SUBROUTINE sbc_cpl_ice_flx( p_frld , palbi   , psst    , pist  )
1600      REAL(wp), INTENT(in   ), DIMENSION(:,:  ) ::   p_frld     ! lead fraction                [0 to 1]
1601      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1602      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature      [Celcius]
1603      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature      [Kelvin]
1604      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) 
1605   END SUBROUTINE sbc_cpl_ice_flx
1606   
1607#endif
1608
1609   !!======================================================================
1610END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.