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

Last change on this file since 2086 was 2086, checked in by omamce, 11 years ago

O.M. : bug correction. Temporary

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