source: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbccpl.F90 @ 2189

Last change on this file since 2189 was 2189, checked in by omamce, 8 years ago

O.M. : small changes to fit IPSLCM5A / NEMO 3.2

File size: 95.9 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: sbccpl.F90 4162 2013-11-07 10:19:49Z cetlod $
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!$$ O.M.         IF ( k_ice <= 1 )  srcv(jpr_ivep)%laction = .FALSE.
397         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
398      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
399      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
400      END SELECT
401
402      !                                                      ! ------------------------- !
403      !                                                      !     Runoffs & Calving     !   
404      !                                                      ! ------------------------- !
405      srcv(jpr_rnf   )%clname = 'O_Runoff'   ;   IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' )   srcv(jpr_rnf)%laction = .TRUE.
406! This isn't right - really just want ln_rnf_emp changed
407!                                                 IF( TRIM( sn_rcv_rnf%cldes ) == 'climato' )   THEN   ;   ln_rnf = .TRUE.
408!                                                 ELSE                                                 ;   ln_rnf = .FALSE.
409!                                                 ENDIF
410      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
411
412      !                                                      ! ------------------------- !
413      !                                                      !    non solar radiation    !   Qns
414      !                                                      ! ------------------------- !
415      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
416      srcv(jpr_qnsice)%clname = 'O_QnsIce'
417      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
418      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
419      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
420      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
421      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
422      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
423      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
424      END SELECT
425      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
426         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
427      !                                                      ! ------------------------- !
428      !                                                      !    solar radiation        !   Qsr
429      !                                                      ! ------------------------- !
430      srcv(jpr_qsroce)%clname = 'O_QsrOce'
431      srcv(jpr_qsrice)%clname = 'O_QsrIce'
432      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
433      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
434      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
435      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
436      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
437      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
438      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
439      END SELECT
440      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
441         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
442      !                                                      ! ------------------------- !
443      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
444      !                                                      ! ------------------------- !
445      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
446      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
447      !
448      ! non solar sensitivity mandatory for LIM ice model
449      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4) &
450         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
451      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
452      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
453         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
454      !                                                      ! ------------------------- !
455      !                                                      !      10m wind module      !   
456      !                                                      ! ------------------------- !
457      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
458      !
459      !                                                      ! ------------------------- !
460      !                                                      !   wind stress module      !   
461      !                                                      ! ------------------------- !
462      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
463      lhftau = srcv(jpr_taum)%laction
464
465      !                                                      ! ------------------------- !
466      !                                                      !      Atmospheric CO2      !
467      !                                                      ! ------------------------- !
468      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
469      !                                                      ! ------------------------- !
470      !                                                      !   topmelt and botmelt     !   
471      !                                                      ! ------------------------- !
472      srcv(jpr_topm )%clname = 'OTopMlt'
473      srcv(jpr_botm )%clname = 'OBotMlt'
474      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
475         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
476            srcv(jpr_topm:jpr_botm)%nct = jpl
477         ELSE
478            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
479         ENDIF
480         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
481      ENDIF
482
483      ! Allocate all parts of frcv used for received fields
484      DO jn = 1, jprcv
485         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
486      END DO
487      ! Allocate taum part of frcv which is used even when not received as coupling field
488      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jn)%nct) )
489      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
490      IF( k_ice /= 0 ) THEN
491         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jn)%nct) )
492         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jn)%nct) )
493      END IF
494
495      ! ================================ !
496      !     Define the send interface    !
497      ! ================================ !
498      ! for each field: define the OASIS name                           (ssnd(:)%clname)
499      !                 define send or not from the namelist parameters (ssnd(:)%laction)
500      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
501     
502      ! default definitions of nsnd
503      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
504         
505      !                                                      ! ------------------------- !
506      !                                                      !    Surface temperature    !
507      !                                                      ! ------------------------- !
508      ssnd(jps_toce)%clname = 'O_SSTSST'
509      ssnd(jps_tice)%clname = 'O_TepIce'
510      ssnd(jps_tmix)%clname = 'O_TepMix'
511      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
512      CASE( 'none'         )       ! nothing to do
513      CASE( 'oce only'             )   ;   ssnd(   jps_toce             )%laction = .TRUE.
514      CASE( 'weighted oce and ice' )
515         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
516         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
517      CASE( 'mixed oce-ice'        )   ;   ssnd(   jps_tmix             )%laction = .TRUE.
518      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
519      END SELECT
520     
521      !                                                      ! ------------------------- !
522      !                                                      !          Albedo           !
523      !                                                      ! ------------------------- !
524      ssnd(jps_albice)%clname = 'O_AlbIce' 
525      ssnd(jps_albmix)%clname = 'O_AlbMix'
526      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
527      CASE( 'none'          )       ! nothing to do
528      CASE( 'weighted ice'  )   ;   ssnd(jps_albice)%laction = .TRUE.
529      CASE( 'mixed oce-ice' )   ;   ssnd(jps_albmix)%laction = .TRUE.
530      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
531      END SELECT
532      !
533      ! Need to calculate oceanic albedo if
534      !     1. sending mixed oce-ice albedo or
535      !     2. receiving mixed oce-ice solar radiation
536      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
537         CALL albedo_oce( zaos, zacs )
538         ! Due to lack of information on nebulosity : mean clear/overcast sky
539         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
540      ENDIF
541
542      !                                                      ! ------------------------- !
543      !                                                      !  Ice fraction & Thickness !
544      !                                                      ! ------------------------- !
545      ssnd(jps_fice)%clname = 'OIceFrc'
546      ssnd(jps_hice)%clname = 'OIceTck'
547      ssnd(jps_hsnw)%clname = 'OSnwTck'
548      IF( k_ice /= 0 ) THEN
549         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
550! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
551         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
552      ENDIF
553
554      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
555      CASE( 'none'         )       ! nothing to do
556      CASE( 'ice and snow' ) 
557         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
558         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
559            ssnd(jps_hice:jps_hsnw)%nct = jpl
560         ELSE
561            IF ( jpl > 1 ) THEN
562CALL ctl_stop( 'sbc_cpl_init: use weighted ice and snow option for sn_snd_thick%cldes if not exchanging category fields' )
563            ENDIF
564         ENDIF
565      CASE ( 'weighted ice and snow' ) 
566         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
567         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
568      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
569      END SELECT
570
571      !                                                      ! ------------------------- !
572      !                                                      !      Surface current      !
573      !                                                      ! ------------------------- !
574      !        ocean currents              !            ice velocities
575      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
576      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
577      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
578      !
579      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
580
581      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
582         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
583      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
584         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
585         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
586      ENDIF
587      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
588      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
589      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
590      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
591      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
592      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
593      CASE( 'weighted oce and ice' )   !   nothing to do
594      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
595      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
596      END SELECT
597
598      !                                                      ! ------------------------- !
599      !                                                      !          CO2 flux         !
600      !                                                      ! ------------------------- !
601      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
602      !
603      ! ================================ !
604      !   initialisation of the coupler  !
605      ! ================================ !
606
607      CALL cpl_prism_define(jprcv, jpsnd)           
608      !
609      IF( ln_dm2dc .AND. ( cpl_prism_freq( jpr_qsroce ) + cpl_prism_freq( jpr_qsrmix ) /= 86400 ) )   &
610         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
611
612      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
613      !
614      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
615      !
616   END SUBROUTINE sbc_cpl_init
617
618
619   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
620      !!----------------------------------------------------------------------
621      !!             ***  ROUTINE sbc_cpl_rcv  ***
622      !!
623      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
624      !!                provide the ocean heat and freshwater fluxes.
625      !!
626      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
627      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
628      !!                to know if the field was really received or not
629      !!
630      !!              --> If ocean stress was really received:
631      !!
632      !!                  - transform the received ocean stress vector from the received
633      !!                 referential and grid into an atmosphere-ocean stress in
634      !!                 the (i,j) ocean referencial and at the ocean velocity point.
635      !!                    The received stress are :
636      !!                     - defined by 3 components (if cartesian coordinate)
637      !!                            or by 2 components (if spherical)
638      !!                     - oriented along geographical   coordinate (if eastward-northward)
639      !!                            or  along the local grid coordinate (if local grid)
640      !!                     - given at U- and V-point, resp.   if received on 2 grids
641      !!                            or at T-point               if received on 1 grid
642      !!                    Therefore and if necessary, they are successively
643      !!                  processed in order to obtain them
644      !!                     first  as  2 components on the sphere
645      !!                     second as  2 components oriented along the local grid
646      !!                     third  as  2 components on the U,V grid
647      !!
648      !!              -->
649      !!
650      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
651      !!             and total ocean freshwater fluxes 
652      !!
653      !! ** Method  :   receive all fields from the atmosphere and transform
654      !!              them into ocean surface boundary condition fields
655      !!
656      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
657      !!                        taum, wndm   wind stres and wind speed module at T-point
658      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
659      !!                                     and the latent heat flux of solid precip. melting
660      !!                        qsr          solar ocean heat fluxes   (ocean only case)
661      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
662      !!----------------------------------------------------------------------
663      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
664      INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation
665      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
666      !!
667      LOGICAL ::    llnewtx, llnewtau      ! update wind stress components and module??
668      INTEGER  ::   ji, jj, jn             ! dummy loop indices
669      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
670      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
671      REAL(wp) ::   zcoef                  ! temporary scalar
672      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
673      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
674      REAL(wp) ::   zzx, zzy               ! temporary variables
675      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
676      !!----------------------------------------------------------------------
677      !
678      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
679      !
680      CALL wrk_alloc( jpi,jpj, ztx, zty )
681
682      IF( kt == nit000 )   CALL sbc_cpl_init( k_ice )          ! initialisation
683
684      !                                                 ! Receive all the atmos. fields (including ice information)
685      isec = ( kt - nit000 ) * NINT( rdttra(1) )             ! date of exchanges
686      DO jn = 1, jprcv                                       ! received fields sent by the atmosphere
687         IF( srcv(jn)%laction )   CALL cpl_prism_rcv( jn, isec, frcv(jn)%z3, nrcvinfo(jn) )
688      END DO
689
690      !                                                      ! ========================= !
691      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
692         !                                                   ! ========================= !
693         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
694         ! => need to be done only when we receive the field
695         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
696            !
697            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
698               !                                                       ! (cartesian to spherical -> 3 to 2 components)
699               !
700               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
701                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
702               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
703               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
704               !
705               IF( srcv(jpr_otx2)%laction ) THEN
706                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
707                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
708                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
709                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
710               ENDIF
711               !
712            ENDIF
713            !
714            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
715               !                                                       ! (geographical to local grid -> rotate the components)
716               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
717               IF( srcv(jpr_otx2)%laction ) THEN
718                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
719               ELSE     
720                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
721               ENDIF
722               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
723               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
724            ENDIF
725            !                             
726            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
727               DO jj = 2, jpjm1                                          ! T ==> (U,V)
728                  DO ji = fs_2, fs_jpim1   ! vector opt.
729                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
730                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
731                  END DO
732               END DO
733               CALL lbc_lnk( frcv(jpr_otx1)%z3(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
734            ENDIF
735            llnewtx = .TRUE.
736         ELSE
737            llnewtx = .FALSE.
738         ENDIF
739         !                                                   ! ========================= !
740      ELSE                                                   !   No dynamical coupling   !
741         !                                                   ! ========================= !
742         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
743         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
744         llnewtx = .TRUE.
745         !
746      ENDIF
747     
748      !                                                      ! ========================= !
749      !                                                      !    wind stress module     !   (taum)
750      !                                                      ! ========================= !
751      !
752      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
753         ! => need to be done only when otx1 was changed
754         IF( llnewtx ) THEN
755!CDIR NOVERRCHK
756            DO jj = 2, jpjm1
757!CDIR NOVERRCHK
758               DO ji = fs_2, fs_jpim1   ! vect. opt.
759                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
760                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
761                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
762               END DO
763            END DO
764            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
765            llnewtau = .TRUE.
766         ELSE
767            llnewtau = .FALSE.
768         ENDIF
769      ELSE
770         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
771         ! Stress module can be negative when received (interpolation problem)
772         IF( llnewtau ) THEN
773            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
774         ENDIF
775      ENDIF
776     
777      !                                                      ! ========================= !
778      !                                                      !      10 m wind speed      !   (wndm)
779      !                                                      ! ========================= !
780      !
781      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
782         ! => need to be done only when taumod was changed
783         IF( llnewtau ) THEN
784            zcoef = 1. / ( zrhoa * zcdrag ) 
785!CDIR NOVERRCHK
786            DO jj = 1, jpj
787!CDIR NOVERRCHK
788               DO ji = 1, jpi 
789                  wndm(ji,jj) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
790               END DO
791            END DO
792         ENDIF
793      ELSE
794         IF ( nrcvinfo(jpr_w10m) == OASIS_Rcv ) wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
795      ENDIF
796
797      ! u(v)tau and taum will be modified by ice model
798      ! -> need to be reset before each call of the ice/fsbc     
799      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
800         !
801         utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
802         vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
803         taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
804         CALL iom_put( "taum_oce", taum )   ! output wind stress module
805         
806      ENDIF
807
808#if defined key_cpl_carbon_cycle
809      !                                                              ! atmosph. CO2 (ppm)
810      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
811#endif
812
813      !                                                      ! ========================= !
814      IF( k_ice <= 1 ) THEN                                  !  heat & freshwater fluxes ! (Ocean only case)
815         !                                                   ! ========================= !
816         !
817         !                                                       ! total freshwater fluxes over the ocean (emp)
818         SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
819         CASE( 'conservative' )
820            emp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
821         CASE( 'oce only', 'oce and ice' )
822            emp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
823         CASE default
824            CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
825         END SELECT
826         !
827         !                                                        ! runoffs and calving (added in emp)
828         IF( srcv(jpr_rnf)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
829         IF( srcv(jpr_cal)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_cal)%z3(:,:,1)
830         !
831!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
832!!gm                                       at least should be optional...
833!!         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN     ! add to the total freshwater budget
834!!            ! remove negative runoff
835!!            zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
836!!            zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
837!!            IF( lk_mpp )   CALL mpp_sum( zcumulpos )   ! sum over the global domain
838!!            IF( lk_mpp )   CALL mpp_sum( zcumulneg )
839!!            IF( zcumulpos /= 0. ) THEN                 ! distribute negative runoff on positive runoff grid points
840!!               zcumulneg = 1.e0 + zcumulneg / zcumulpos
841!!               frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
842!!            ENDIF     
843!!            ! add runoff to e-p
844!!            emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
845!!         ENDIF
846!!gm  end of internal cooking
847         !
848         !                                                       ! non solar heat flux over the ocean (qns)
849         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
850         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
851         ! add the latent heat of solid precip. melting
852         IF( srcv(jpr_snow  )%laction )   THEN                         ! update qns over the free ocean with:
853              qns(:,:) = qns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus  & ! energy for melting solid precipitation over the free ocean
854           &           - emp(:,:) * sst_m(:,:) * rcp                   ! remove heat content due to mass flux (assumed to be at SST)
855         ENDIF
856
857         !                                                       ! solar flux over the ocean          (qsr)
858         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
859         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
860         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
861         !
862 
863      ENDIF
864      !
865      CALL wrk_dealloc( jpi,jpj, ztx, zty )
866      !
867      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
868      !
869   END SUBROUTINE sbc_cpl_rcv
870   
871
872   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
873      !!----------------------------------------------------------------------
874      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
875      !!
876      !! ** Purpose :   provide the stress over sea-ice in coupled mode
877      !!
878      !! ** Method  :   transform the received stress from the atmosphere into
879      !!             an atmosphere-ice stress in the (i,j) ocean referencial
880      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
881      !!                'C'-grid : i- (j-) components given at U- (V-) point
882      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
883      !!
884      !!                The received stress are :
885      !!                 - defined by 3 components (if cartesian coordinate)
886      !!                        or by 2 components (if spherical)
887      !!                 - oriented along geographical   coordinate (if eastward-northward)
888      !!                        or  along the local grid coordinate (if local grid)
889      !!                 - given at U- and V-point, resp.   if received on 2 grids
890      !!                        or at a same point (T or I) if received on 1 grid
891      !!                Therefore and if necessary, they are successively
892      !!             processed in order to obtain them
893      !!                 first  as  2 components on the sphere
894      !!                 second as  2 components oriented along the local grid
895      !!                 third  as  2 components on the cp_ice_msh point
896      !!
897      !!                Except in 'oce and ice' case, only one vector stress field
898      !!             is received. It has already been processed in sbc_cpl_rcv
899      !!             so that it is now defined as (i,j) components given at U-
900      !!             and V-points, respectively. Therefore, only the third
901      !!             transformation is done and only if the ice-grid is a 'I'-grid.
902      !!
903      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
904      !!----------------------------------------------------------------------
905      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
906      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
907      !!
908      INTEGER ::   ji, jj                          ! dummy loop indices
909      INTEGER ::   itx                             ! index of taux over ice
910      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
911      !!----------------------------------------------------------------------
912      !
913      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
914      !
915      CALL wrk_alloc( jpi,jpj, ztx, zty )
916
917!AC Pour eviter un stress nul sur la glace dans le cas mixed oce-ice
918      IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN   ;   itx =  jpr_itx1   
919      ELSE                                ;   itx =  jpr_otx1
920      ENDIF
921
922      ! do something only if we just received the stress from atmosphere
923      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
924
925         !                                                                                              ! ======================= !
926!AC Pour eviter un stress nul sur la glace dans le cas mixes oce-ice
927         IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN               !   ice stress received   !
928            !                                                                                           ! ======================= !
929           
930            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
931               !                                                       ! (cartesian to spherical -> 3 to 2 components)
932               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
933                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
934               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
935               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
936               !
937               IF( srcv(jpr_itx2)%laction ) THEN
938                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
939                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
940                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
941                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
942               ENDIF
943               !
944            ENDIF
945            !
946            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
947               !                                                       ! (geographical to local grid -> rotate the components)
948               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
949               IF( srcv(jpr_itx2)%laction ) THEN
950                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
951               ELSE
952                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
953               ENDIF
954               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
955               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
956            ENDIF
957            !                                                   ! ======================= !
958         ELSE                                                   !     use ocean stress    !
959            !                                                   ! ======================= !
960            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
961            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
962            !
963         ENDIF
964
965         !                                                      ! ======================= !
966         !                                                      !     put on ice grid     !
967         !                                                      ! ======================= !
968         !   
969         !                                                  j+1   j     -----V---F
970         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
971         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
972         !                                                               |       |
973         !                                                   j    j-1   -I-------|
974         !                                               (for I)         |       |
975         !                                                              i-1  i   i
976         !                                                               i      i+1 (for I)
977         SELECT CASE ( cp_ice_msh )
978            !
979         CASE( 'I' )                                         ! B-grid ==> I
980            SELECT CASE ( srcv(jpr_itx1)%clgrid )
981            CASE( 'U' )
982               DO jj = 2, jpjm1                                   ! (U,V) ==> I
983                  DO ji = 2, jpim1   ! NO vector opt.
984                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
985                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
986                  END DO
987               END DO
988            CASE( 'F' )
989               DO jj = 2, jpjm1                                   ! F ==> I
990                  DO ji = 2, jpim1   ! NO vector opt.
991                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
992                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
993                  END DO
994               END DO
995            CASE( 'T' )
996               DO jj = 2, jpjm1                                   ! T ==> I
997                  DO ji = 2, jpim1   ! NO vector opt.
998                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
999                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
1000                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
1001                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1002                  END DO
1003               END DO
1004            CASE( 'I' )
1005               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
1006               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1007            END SELECT
1008            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1009               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1010            ENDIF
1011            !
1012         CASE( 'F' )                                         ! B-grid ==> F
1013            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1014            CASE( 'U' )
1015               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1016                  DO ji = fs_2, fs_jpim1   ! vector opt.
1017                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1018                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1019                  END DO
1020               END DO
1021            CASE( 'I' )
1022               DO jj = 2, jpjm1                                   ! I ==> F
1023                  DO ji = 2, jpim1   ! NO vector opt.
1024                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1025                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1026                  END DO
1027               END DO
1028            CASE( 'T' )
1029               DO jj = 2, jpjm1                                   ! T ==> F
1030                  DO ji = 2, jpim1   ! NO vector opt.
1031                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1032                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1033                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1034                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1035                  END DO
1036               END DO
1037            CASE( 'F' )
1038               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1039               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1040            END SELECT
1041            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1042               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1043            ENDIF
1044            !
1045         CASE( 'C' )                                         ! C-grid ==> U,V
1046            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1047            CASE( 'U' )
1048               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1049               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1050            CASE( 'F' )
1051               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1052                  DO ji = fs_2, fs_jpim1   ! vector opt.
1053                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1054                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1055                  END DO
1056               END DO
1057            CASE( 'T' )
1058               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1059                  DO ji = fs_2, fs_jpim1   ! vector opt.
1060                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1061                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1062                  END DO
1063               END DO
1064            CASE( 'I' )
1065               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1066                  DO ji = 2, jpim1   ! NO vector opt.
1067                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1068                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1069                  END DO
1070               END DO
1071            END SELECT
1072            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1073               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1074            ENDIF
1075         END SELECT
1076
1077      ENDIF
1078      !   
1079      CALL wrk_dealloc( jpi,jpj, ztx, zty )
1080      !
1081      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
1082      !
1083   END SUBROUTINE sbc_cpl_ice_tau
1084   
1085
1086   SUBROUTINE sbc_cpl_ice_flx( p_frld  , palbi   , psst    , pist    )
1087      !!----------------------------------------------------------------------
1088      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1089      !!
1090      !! ** Purpose :   provide the heat and freshwater fluxes of the
1091      !!              ocean-ice system.
1092      !!
1093      !! ** Method  :   transform the fields received from the atmosphere into
1094      !!             surface heat and fresh water boundary condition for the
1095      !!             ice-ocean system. The following fields are provided:
1096      !!              * total non solar, solar and freshwater fluxes (qns_tot,
1097      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1098      !!             NB: emp_tot include runoffs and calving.
1099      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1100      !!             emp_ice = sublimation - solid precipitation as liquid
1101      !!             precipitation are re-routed directly to the ocean and
1102      !!             runoffs and calving directly enter the ocean.
1103      !!              * solid precipitation (sprecip), used to add to qns_tot
1104      !!             the heat lost associated to melting solid precipitation
1105      !!             over the ocean fraction.
1106      !!       ===>> CAUTION here this changes the net heat flux received from
1107      !!             the atmosphere
1108      !!
1109      !!                  - the fluxes have been separated from the stress as
1110      !!                 (a) they are updated at each ice time step compare to
1111      !!                 an update at each coupled time step for the stress, and
1112      !!                 (b) the conservative computation of the fluxes over the
1113      !!                 sea-ice area requires the knowledge of the ice fraction
1114      !!                 after the ice advection and before the ice thermodynamics,
1115      !!                 so that the stress is updated before the ice dynamics
1116      !!                 while the fluxes are updated after it.
1117      !!
1118      !! ** Action  :   update at each nf_ice time step:
1119      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1120      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1121      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
1122      !!                   emp_ice            ice sublimation - solid precipitation over the ice
1123      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
1124      !!                   sprecip             solid precipitation over the ocean 
1125      !!----------------------------------------------------------------------
1126      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1127      ! optional arguments, used only in 'mixed oce-ice' case
1128      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1129      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celcius]
1130      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature     [Kelvin]
1131      !
1132      INTEGER ::   jl   ! dummy loop index
1133      REAL(wp), POINTER, DIMENSION(:,:) ::   zcptn, ztmp, zicefr
1134      !!----------------------------------------------------------------------
1135      !
1136      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1137      !
1138      CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr )
1139
1140      zicefr(:,:) = 1.- p_frld(:,:)
1141      zcptn(:,:) = rcp * sst_m(:,:)
1142      !
1143      !                                                      ! ========================= !
1144      !                                                      !    freshwater budget      !   (emp)
1145      !                                                      ! ========================= !
1146      !
1147      !                                                           ! total Precipitations - total Evaporation (emp_tot)
1148      !                                                           ! solid precipitation  - sublimation       (emp_ice)
1149      !                                                           ! solid Precipitation                      (sprecip)
1150      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1151      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1152         sprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                 ! May need to ensure positive here
1153         tprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + sprecip (:,:) ! May need to ensure positive here
1154         emp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - tprecip(:,:)
1155         emp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1156                           CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation
1157         IF( lk_diaar5 )   CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
1158         ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1159                           CALL iom_put( 'evap_ao_cea'  , ztmp                            )   ! ice-free oce evap (cell average)
1160         IF( lk_diaar5 )   CALL iom_put( 'hflx_evap_cea', ztmp(:,:         ) * zcptn(:,:) )   ! heat flux from from evap (cell ave)
1161      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1162         emp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1163         emp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
1164         sprecip(:,:) = - frcv(jpr_semp)%z3(:,:,1) + frcv(jpr_ievp)%z3(:,:,1)
1165      END SELECT
1166
1167      CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
1168      CALL iom_put( 'snow_ao_cea', sprecip(:,:         ) * p_frld(:,:)    )   ! Snow        over ice-free ocean  (cell average)
1169      CALL iom_put( 'snow_ai_cea', sprecip(:,:         ) * zicefr(:,:)    )   ! Snow        over sea-ice         (cell average)
1170      CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
1171      !   
1172      !                                                           ! runoffs and calving (put in emp_tot)
1173      IF( srcv(jpr_rnf)%laction ) THEN
1174         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)
1175                           CALL iom_put( 'runoffs'      , frcv(jpr_rnf)%z3(:,:,1)              )   ! rivers
1176         IF( lk_diaar5 )   CALL iom_put( 'hflx_rnf_cea' , frcv(jpr_rnf)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from rivers
1177      ENDIF
1178      IF( srcv(jpr_cal)%laction ) THEN
1179         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1180         CALL iom_put( 'calving', frcv(jpr_cal)%z3(:,:,1) )
1181      ENDIF
1182      !
1183!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
1184!!gm                                       at least should be optional...
1185!!       ! remove negative runoff                            ! sum over the global domain
1186!!       zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
1187!!       zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
1188!!       IF( lk_mpp )   CALL mpp_sum( zcumulpos )
1189!!       IF( lk_mpp )   CALL mpp_sum( zcumulneg )
1190!!       IF( zcumulpos /= 0. ) THEN                          ! distribute negative runoff on positive runoff grid points
1191!!          zcumulneg = 1.e0 + zcumulneg / zcumulpos
1192!!          frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
1193!!       ENDIF     
1194!!       emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)   ! add runoff to e-p
1195!!
1196!!gm  end of internal cooking
1197
1198      !                                                      ! ========================= !
1199      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1200      !                                                      ! ========================= !
1201      CASE( 'oce only' )                                     ! the required field is directly provided
1202         qns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
1203      CASE( 'conservative' )                                      ! the required fields are directly provided
1204         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1205         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1206            qns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1207         ELSE
1208            ! Set all category values equal for the moment
1209            DO jl=1,jpl
1210               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1211            ENDDO
1212         ENDIF
1213      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
1214         qns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1215         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1216            DO jl=1,jpl
1217               qns_tot(:,:   ) = qns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1218               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1219            ENDDO
1220         ELSE
1221            DO jl=1,jpl
1222               qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1223               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1224            ENDDO
1225         ENDIF
1226      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
1227! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1228         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1229         qns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1230            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1231            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
1232      END SELECT
1233      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus
1234      qns_tot(:,:) = qns_tot(:,:)                         &            ! qns_tot update over free ocean with:
1235         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
1236         &          - (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
1237         &             - emp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
1238      IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1239!!gm
1240!!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in
1241!!    the flux that enter the ocean....
1242!!    moreover 1 - it is not diagnose anywhere....
1243!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1244!!
1245!! similar job should be done for snow and precipitation temperature
1246      !                                     
1247      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting
1248         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting
1249         qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)
1250         IF( lk_diaar5 )   CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1251      ENDIF
1252
1253      !                                                      ! ========================= !
1254      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1255      !                                                      ! ========================= !
1256      CASE( 'oce only' )
1257         qsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
1258      CASE( 'conservative' )
1259         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1260         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1261            qsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1262         ELSE
1263            ! Set all category values equal for the moment
1264            DO jl=1,jpl
1265               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1266            ENDDO
1267         ENDIF
1268         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1269         qsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1270      CASE( 'oce and ice' )
1271         qsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1272         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1273            DO jl=1,jpl
1274               qsr_tot(:,:   ) = qsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1275               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1276            ENDDO
1277         ELSE
1278            DO jl=1,jpl
1279               qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1280               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1281            ENDDO
1282         ENDIF
1283      CASE( 'mixed oce-ice' )
1284         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1285! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1286!       Create solar heat flux over ice using incoming solar heat flux and albedos
1287!       ( see OASIS3 user guide, 5th edition, p39 )
1288         qsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1289            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1290            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1291      END SELECT
1292      IF( ln_dm2dc ) THEN   ! modify qsr to include the diurnal cycle
1293         qsr_tot(:,:  ) = sbc_dcy( qsr_tot(:,:  ) )
1294         DO jl=1,jpl
1295            qsr_ice(:,:,jl) = sbc_dcy( qsr_ice(:,:,jl) )
1296         ENDDO
1297      ENDIF
1298
1299      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )
1300      CASE ('coupled')
1301         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1302            dqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
1303         ELSE
1304            ! Set all category values equal for the moment
1305            DO jl=1,jpl
1306               dqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
1307            ENDDO
1308         ENDIF
1309      END SELECT
1310
1311      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )
1312      CASE ('coupled')
1313         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
1314         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
1315      END SELECT
1316
1317      !    Ice Qsr penetration used (only?)in lim2 or lim3
1318      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer
1319      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
1320      ! Coupled case: since cloud cover is not received from atmosphere
1321      !               ===> defined as constant value -> definition done in sbc_cpl_init
1322      fr1_i0(:,:) = 0.18
1323      fr2_i0(:,:) = 0.82
1324
1325
1326      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
1327      !
1328      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
1329      !
1330   END SUBROUTINE sbc_cpl_ice_flx
1331   
1332   
1333   SUBROUTINE sbc_cpl_snd( kt )
1334      !!----------------------------------------------------------------------
1335      !!             ***  ROUTINE sbc_cpl_snd  ***
1336      !!
1337      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
1338      !!
1339      !! ** Method  :   send to the atmosphere through a call to cpl_prism_snd
1340      !!              all the needed fields (as defined in sbc_cpl_init)
1341      !!----------------------------------------------------------------------
1342      INTEGER, INTENT(in) ::   kt
1343      !
1344      INTEGER ::   ji, jj, jl   ! dummy loop indices
1345      INTEGER ::   isec, info   ! local integer
1346      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
1347      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
1348      !!----------------------------------------------------------------------
1349      !
1350      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
1351      !
1352      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1353      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1354
1355      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
1356
1357      zfr_l(:,:) = 1.- fr_i(:,:)
1358
1359      !                                                      ! ------------------------- !
1360      !                                                      !    Surface temperature    !   in Kelvin
1361      !                                                      ! ------------------------- !
1362      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
1363         SELECT CASE( sn_snd_temp%cldes)
1364         CASE( 'oce only'             )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0
1365         CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:)   
1366            SELECT CASE( sn_snd_temp%clcat )
1367            CASE( 'yes' )   
1368               ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1369            CASE( 'no' )
1370               ztmp3(:,:,:) = 0.0
1371               DO jl=1,jpl
1372                  ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
1373               ENDDO
1374            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1375            END SELECT
1376         CASE( 'mixed oce-ice'        )   
1377            ztmp1(:,:) = ( tsn(:,:,1,1) + rt0 ) * zfr_l(:,:) 
1378            DO jl=1,jpl
1379               ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
1380            ENDDO
1381         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
1382         END SELECT
1383         IF( ssnd(jps_toce)%laction )   CALL cpl_prism_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1384         IF( ssnd(jps_tice)%laction )   CALL cpl_prism_snd( jps_tice, isec, ztmp3, info )
1385         IF( ssnd(jps_tmix)%laction )   CALL cpl_prism_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1386      ENDIF
1387      !
1388      !                                                      ! ------------------------- !
1389      !                                                      !           Albedo          !
1390      !                                                      ! ------------------------- !
1391      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
1392         ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1393         CALL cpl_prism_snd( jps_albice, isec, ztmp3, info )
1394      ENDIF
1395      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
1396         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
1397         DO jl=1,jpl
1398            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
1399         ENDDO
1400         CALL cpl_prism_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1401      ENDIF
1402      !                                                      ! ------------------------- !
1403      !                                                      !  Ice fraction & Thickness !
1404      !                                                      ! ------------------------- !
1405      ! Send ice fraction field
1406      IF( ssnd(jps_fice)%laction ) THEN
1407         SELECT CASE( sn_snd_thick%clcat )
1408         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
1409         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
1410         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1411         END SELECT
1412         CALL cpl_prism_snd( jps_fice, isec, ztmp3, info )
1413      ENDIF
1414
1415      ! Send ice and snow thickness field
1416      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
1417         SELECT CASE( sn_snd_thick%cldes)
1418         CASE( 'none'                  )       ! nothing to do
1419         CASE( 'weighted ice and snow' )   
1420            SELECT CASE( sn_snd_thick%clcat )
1421            CASE( 'yes' )   
1422               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
1423               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
1424            CASE( 'no' )
1425               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
1426               DO jl=1,jpl
1427                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
1428                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
1429               ENDDO
1430            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1431            END SELECT
1432         CASE( 'ice and snow'         )   
1433            ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
1434            ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
1435         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
1436         END SELECT
1437         IF( ssnd(jps_hice)%laction )   CALL cpl_prism_snd( jps_hice, isec, ztmp3, info )
1438         IF( ssnd(jps_hsnw)%laction )   CALL cpl_prism_snd( jps_hsnw, isec, ztmp4, info )
1439      ENDIF
1440      !
1441#if defined key_cpl_carbon_cycle
1442      !                                                      ! ------------------------- !
1443      !                                                      !  CO2 flux from PISCES     !
1444      !                                                      ! ------------------------- !
1445      IF( ssnd(jps_co2)%laction )   CALL cpl_prism_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
1446      !
1447#endif
1448      !                                                      ! ------------------------- !
1449      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
1450         !                                                   ! ------------------------- !
1451         !   
1452         !                                                  j+1   j     -----V---F
1453         ! surface velocity always sent from T point                     !       |
1454         !                                                        j      |   T   U
1455         !                                                               |       |
1456         !                                                   j    j-1   -I-------|
1457         !                                               (for I)         |       |
1458         !                                                              i-1  i   i
1459         !                                                               i      i+1 (for I)
1460         SELECT CASE( TRIM( sn_snd_crt%cldes ) )
1461         CASE( 'oce only'             )      ! C-grid ==> T
1462            DO jj = 2, jpjm1
1463               DO ji = fs_2, fs_jpim1   ! vector opt.
1464                  zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
1465                  zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
1466               END DO
1467            END DO
1468         CASE( 'weighted oce and ice' )   
1469            SELECT CASE ( cp_ice_msh )
1470            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1471               DO jj = 2, jpjm1
1472                  DO ji = fs_2, fs_jpim1   ! vector opt.
1473                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1474                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
1475                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1476                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1477                  END DO
1478               END DO
1479            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1480               DO jj = 2, jpjm1
1481                  DO ji = 2, jpim1   ! NO vector opt.
1482                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1483                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1484                     zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1485                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1486                     zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1487                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1488                  END DO
1489               END DO
1490            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1491               DO jj = 2, jpjm1
1492                  DO ji = 2, jpim1   ! NO vector opt.
1493                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1494                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1495                     zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1496                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1497                     zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1498                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1499                  END DO
1500               END DO
1501            END SELECT
1502            CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
1503         CASE( 'mixed oce-ice'        )
1504            SELECT CASE ( cp_ice_msh )
1505            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1506               DO jj = 2, jpjm1
1507                  DO ji = fs_2, fs_jpim1   ! vector opt.
1508                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
1509                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1510                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
1511                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1512                  END DO
1513               END DO
1514            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1515               DO jj = 2, jpjm1
1516                  DO ji = 2, jpim1   ! NO vector opt.
1517                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1518                        &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1519                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1520                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1521                        &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1522                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1523                  END DO
1524               END DO
1525            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1526               DO jj = 2, jpjm1
1527                  DO ji = 2, jpim1   ! NO vector opt.
1528                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1529                        &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1530                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1531                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1532                        &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1533                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1534                  END DO
1535               END DO
1536            END SELECT
1537         END SELECT
1538         CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
1539         !
1540         !
1541         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
1542            !                                                                     ! Ocean component
1543            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
1544            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
1545            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
1546            zoty1(:,:) = ztmp2(:,:)
1547            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
1548               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
1549               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
1550               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
1551               zity1(:,:) = ztmp2(:,:)
1552            ENDIF
1553         ENDIF
1554         !
1555         ! spherical coordinates to cartesian -> 2 components to 3 components
1556         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
1557            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
1558            ztmp2(:,:) = zoty1(:,:)
1559            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
1560            !
1561            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
1562               ztmp1(:,:) = zitx1(:,:)
1563               ztmp1(:,:) = zity1(:,:)
1564               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
1565            ENDIF
1566         ENDIF
1567         !
1568         IF( ssnd(jps_ocx1)%laction )   CALL cpl_prism_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
1569         IF( ssnd(jps_ocy1)%laction )   CALL cpl_prism_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
1570         IF( ssnd(jps_ocz1)%laction )   CALL cpl_prism_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
1571         !
1572         IF( ssnd(jps_ivx1)%laction )   CALL cpl_prism_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
1573         IF( ssnd(jps_ivy1)%laction )   CALL cpl_prism_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
1574         IF( ssnd(jps_ivz1)%laction )   CALL cpl_prism_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
1575         !
1576      ENDIF
1577      !
1578      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1579      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1580      !
1581      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
1582      !
1583   END SUBROUTINE sbc_cpl_snd
1584   
1585#else
1586   !!----------------------------------------------------------------------
1587   !!   Dummy module                                            NO coupling
1588   !!----------------------------------------------------------------------
1589   USE par_kind        ! kind definition
1590CONTAINS
1591   SUBROUTINE sbc_cpl_snd( kt )
1592      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt
1593   END SUBROUTINE sbc_cpl_snd
1594   !
1595   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
1596      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt, k_fsbc, k_ice
1597   END SUBROUTINE sbc_cpl_rcv
1598   !
1599   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1600      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1601      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1602      p_taui(:,:) = 0.   ;   p_tauj(:,:) = 0. ! stupid definition to avoid warning message when compiling...
1603      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?'
1604   END SUBROUTINE sbc_cpl_ice_tau
1605   !
1606   SUBROUTINE sbc_cpl_ice_flx( p_frld , palbi   , psst    , pist  )
1607      REAL(wp), INTENT(in   ), DIMENSION(:,:  ) ::   p_frld     ! lead fraction                [0 to 1]
1608      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1609      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature      [Celcius]
1610      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature      [Kelvin]
1611      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) 
1612   END SUBROUTINE sbc_cpl_ice_flx
1613   
1614#endif
1615
1616   !!======================================================================
1617END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.