New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbccpl.F90 in branches/2011/dev_r2802_UKMO8_sbccpl/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2011/dev_r2802_UKMO8_sbccpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 2885

Last change on this file since 2885 was 2885, checked in by charris, 13 years ago

#662 Changes to complement those in [2884] to allow NEMO-CICE model to compile OK.

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