source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 4391

Last change on this file since 4391 was 4391, checked in by charris, 7 years ago

#1206 Fix for bug introduced under #1107.

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