source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 3488

Last change on this file since 3488 was 3488, checked in by acc, 9 years ago

Branch: dev_r3385_NOCS04_HAMF; #665. Stage 3 of 2012 development: Rationalisation of code. Added LIM3 changes, corrected coupled changes and highlighted areas of concern in CICE interface

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