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

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

#662 Fix to ordering of ice and snow thickness fields.

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