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 @ 2832

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

#662 Latest code changes for sbccpl. These include bug-fixes, control of topmelt and botmelt coupling fields for running CICE with the UM, some more steps towards being able to run with multiple categories for LIM3 and a change in the arguments passed to sbc_cpl_ice_flx.

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