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

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

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

Last change on this file since 3039 was 3039, checked in by charris, 12 years ago

#662 Changes to header information + more namelists updated.

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