source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 2708

Last change on this file since 2708 was 2708, checked in by smasson, 10 years ago

dynamic memory: correction for coupled model

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