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.
ocesbc.F90 in trunk/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMO/OPA_SRC/SBC/ocesbc.F90 @ 699

Last change on this file since 699 was 699, checked in by smasson, 17 years ago

insert revision Id

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 36.0 KB
Line 
1MODULE ocesbc
2   !!======================================================================
3   !!                     ***  MODULE  ocesbc  ***
4   !!                     Ocean surface boundary conditions
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   oce_sbc     : ???
9   !!   oce_sbc_dmp : ???
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce            ! dynamics and tracers variables
13   USE dom_oce        ! ocean space domain variables
14   USE cpl_oce        ! coupled ocean-atmosphere variables
15   USE ice_oce        ! sea-ice variable
16   USE blk_oce        ! bulk variables
17   USE flx_oce        ! sea-ice/ocean forcings variables
18   USE phycst         ! Define parameters for the routines
19   USE taumod         ! surface stress forcing
20   USE flxmod         ! thermohaline fluxes
21   USE flxrnf         ! runoffs forcing
22   USE tradmp         ! damping salinity trend
23   USE dtatem         ! ocean temperature data
24   USE dtasal         ! ocean salinity data
25   USE ocfzpt         ! surface ocean freezing point
26   USE lbclnk         ! ocean lateral boundary condition
27   USE lib_mpp        ! distribued memory computing library
28   USE in_out_manager ! I/O manager
29   USE prtctl         ! Print control
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Accessibility
35   PUBLIC oce_sbc    ! routine called by step
36
37   !! * Shared module variables
38   REAL(wp), PUBLIC ::   &  !:
39      aplus, aminus,     &  !:
40      empold = 0.e0         !: current year freshwater budget correction
41   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
42      qt  ,         &  !: total surface heat flux (w/m2)
43      qsr ,         &  !: solar radiation (w/m2)
44      qla ,         &  !: Latent heat flux (W/m2)
45      qlw ,         &  !: Long Wave Heat Flux (W/m2)
46      qsb ,         &  !: Sensible Heat Flux (W/m2)
47      emp ,         &  !: evaporation minus precipitation (kg/m2/s = mm/s)
48      emps,         &  !: evaporation - precipitation (free surface)
49      qrp ,         &  !: heat flux damping (w/m2)
50      erp              !: evaporation damping (kg/m2/s = mm/s)
51#if ! defined key_dynspg_rl
52   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
53      dmp              !: internal dampind term
54#endif
55
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58
59   !!----------------------------------------------------------------------
60   !!  OPA 9.0 , LOCEAN-IPSL (2005)
61   !! $Id$
62   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
63   !!----------------------------------------------------------------------
64CONTAINS
65
66#if defined key_ice_lim
67   !!----------------------------------------------------------------------
68   !!   'key_ice_lim' :                                   LIM sea-ice model
69   !!----------------------------------------------------------------------
70# if defined key_coupled
71      !!----------------------------------------------------------------------
72      !!   'key_coupled' :                            Coupled Ocean/Atmosphere
73      !!----------------------------------------------------------------------
74
75   SUBROUTINE oce_sbc( kt )
76      !!---------------------------------------------------------------------
77      !!                 ***  ROUTINE oce_sbc  ***
78      !!                   
79      !! ** Purpose :   Ocean surface boundaries conditions with
80      !!        Louvain la Neuve Sea Ice Model in coupled mode
81      !!
82      !! History :
83      !!   1.0  !  00-10  (O. Marti)  Original code
84      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
85      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
86      !!----------------------------------------------------------------------
87      !! * Arguments
88      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
89
90      !! * Local declarations
91      INTEGER ::   ji, jj                   ! dummy loop indices
92      REAL(wp) ::   ztx, ztaux, zty, ztauy
93      REAL(wp) ::   ztdta, ztgel, zqrp
94      !!----------------------------------------------------------------------
95
96      ! 1. initialization to zero at kt = nit000
97      ! ---------------------------------------
98
99      IF( kt == nit000 ) THEN     
100         qsr   (:,:) = 0.e0
101         freeze(:,:) = 0.e0
102         qt    (:,:) = 0.e0
103         qrp   (:,:) = 0.e0
104         emp   (:,:) = 0.e0
105         emps  (:,:) = 0.e0
106         erp   (:,:) = 0.e0
107#if ! defined key_dynspg_rl
108         dmp   (:,:) = 0.e0
109#endif
110      ENDIF
111
112      IF( MOD( kt-1, nfice ) == 0 ) THEN
113
114         CALL oce_sbc_dmp   ! Computation of internal and evaporation damping terms       
115
116         ! Surface heat flux (W/m2)
117         ! -----------------------
118
119         ! restoring heat flux
120         DO jj = 1, jpj
121            DO ji = 1, jpi
122               ztgel = fzptn(ji,jj)
123#if defined key_dtasst
124               ztdta = MAX( sst(ji,jj),    ztgel )
125#else
126               ztdta = MAX( t_dta(ji,jj,1), ztgel )
127#endif
128               zqrp = dqdt0 * ( tb(ji,jj,1) - ztdta )
129
130               qrp(ji,jj) = (1.0-freeze(ji,jj) ) * zqrp
131            END DO
132         END DO
133
134         ! non solar heat flux + solar flux + restoring
135         qt (:,:) = fnsolar(:,:) + fsolar(:,:) + qrp(:,:)
136
137         ! solar flux
138         qsr(:,:) = fsolar(:,:)
139
140#if ! defined key_dynspg_rl 
141         ! total concentration/dilution effect (use on SSS)
142         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
143
144         ! total volume flux (use on sea-surface height)
145         emp (:,:) = fmass(:,:)  -  dmp(:,:) + runoff(:,:) + erp(:,:)
146#else
147         ! Rigid-lid (emp=emps=E-P-R+Erp)
148         ! freshwater flux
149         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
150         emp (:,:) = emps(:,:)
151#endif
152
153         DO jj = 1, jpjm1
154            DO ji = 1, fs_jpim1   ! vertor opt.
155               ztx   = 0.5 * ( freeze(ji+1,jj) + freeze(ji+1,jj+1) )
156               ztaux = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) )
157               taux(ji,jj) = (1.0-ztx) * taux(ji,jj) + ztx * ztaux
158
159               zty   = 0.5 * ( freeze(ji,jj+1) + freeze(ji+1,jj+1) )
160               ztauy = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) )
161               tauy(ji,jj) = (1.0-zty) * tauy(ji,jj) + zty * ztauy
162            END DO
163         END DO
164         CALL lbc_lnk( taux, 'U', -1. )
165         CALL lbc_lnk( tauy, 'V', -1. )   
166
167         ! Re-initialization of fluxes
168         sst_io(:,:) = 0.e0
169         sss_io(:,:) = 0.e0
170         u_io  (:,:) = 0.e0
171         v_io  (:,:) = 0.e0
172         gtaux (:,:) = 0.e0
173         gtauy (:,:) = 0.e0
174
175      ENDIF
176
177   END SUBROUTINE oce_sbc
178
179# elif defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core
180      !!----------------------------------------------------------------------
181      !!   'key_ice_lim'                              with  LIM sea-ice model
182      !!----------------------------------------------------------------------
183
184   SUBROUTINE oce_sbc( kt )
185      !!---------------------------------------------------------------------
186      !!                   ***  ROUTINE oce_sbc  ***
187      !!                   
188      !! ** Purpose : - Ocean surface boundary conditions with LIM sea-ice
189      !!        model in forced mode using bulk formulea
190      !!
191      !! History :
192      !!   1.0  !  99-11  (M. Imbard)  Original code
193      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
194      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
195      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
196      !!----------------------------------------------------------------------
197      !! * arguments
198      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
199
200      !! * Local declarations
201      INTEGER  ::   ji, jj                   ! dummy loop indices
202      REAL(wp) ::   ztx, ztaux, zty, ztauy
203      !!----------------------------------------------------------------------
204
205      ! 1. initialization to zero at kt = nit000
206      ! ---------------------------------------
207
208      IF( kt == nit000 ) THEN     
209         qsr    (:,:) = 0.e0
210         qt     (:,:) = 0.e0
211         qla    (:,:) = 0.e0
212         qlw    (:,:) = 0.e0
213         qsb    (:,:) = 0.e0
214         qrp    (:,:) = 0.e0
215         emp    (:,:) = 0.e0
216         emps   (:,:) = 0.e0
217         erp    (:,:) = 0.e0
218#if ! defined key_dynspg_rl
219         dmp    (:,:) = 0.e0
220#endif
221      ENDIF
222#if defined key_flx_core
223      qla(:,:) = (1 - freeze(:,:))*qla_oce + freeze(:,:)*qla_ice(:,:)
224      qsb(:,:) = (1 - freeze(:,:))*qsb_oce + freeze(:,:)*qsb_ice(:,:)
225      qlw(:,:) = (1 - freeze(:,:))*qlw_oce + freeze(:,:)*qlw_ice(:,:)
226#endif
227      IF( MOD( kt-1, nfice ) == 0 ) THEN
228
229         CALL oce_sbc_dmp       ! Computation of internal and evaporation damping terms       
230
231         ! Surface Ocean fluxes
232         ! ====================
233
234         ! Surface heat flux (W/m2)
235         ! -----------------
236
237         qt  (:,:) = fnsolar(:,:) + fsolar(:,:)     ! non solar heat flux + solar flux
238         qsr (:,:) = fsolar(:,:)                     ! solar flux
239
240#if ! defined key_dynspg_rl     
241         ! total concentration/dilution effect (use on SSS)
242         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) + empold
243
244         ! total volume flux (use on sea-surface height)
245         emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold     
246#else
247         ! Rigid-lid (emp=emps=E-P-R+Erp)
248         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)     ! freshwater flux
249         emp (:,:) = emps(:,:)
250
251#endif
252
253         ! Surface stress
254         ! --------------
255
256         ! update the stress beloww sea-ice area
257         DO jj = 1, jpjm1
258            DO ji = 1, fs_jpim1   ! vertor opt.
259               ztx         = MAX( freezn(ji,jj), freezn(ji,jj+1) )   ! ice/ocean indicator at U- and V-points
260               zty         = MAX( freezn(ji,jj), freezn(ji+1,jj) )
261               ztaux       = 0.5 *( ftaux(ji+1,jj) + ftaux(ji+1,jj+1) ) ! ice-ocean stress at U- and V-points
262               ztauy       = 0.5 *( ftauy(ji,jj+1) + ftauy(ji+1,jj+1) )
263               taux(ji,jj) = (1.-ztx) * taux(ji,jj) + ztx * ztaux    ! stress at the ocean surface
264               tauy(ji,jj) = (1.-zty) * tauy(ji,jj) + zty * ztauy
265            END DO
266         END DO
267
268         ! boundary condition on the stress (taux,tauy)
269         CALL lbc_lnk( taux, 'U', -1. )
270         CALL lbc_lnk( tauy, 'V', -1. )
271
272         ! Re-initialization of fluxes
273         sst_io(:,:) = 0.e0
274         sss_io(:,:) = 0.e0
275         u_io  (:,:) = 0.e0
276         v_io  (:,:) = 0.e0
277
278      ENDIF
279
280   END SUBROUTINE oce_sbc
281
282# else
283      !!----------------------------------------------------------------------
284      !!   Error option               LIM sea-ice model requires bulk formulea
285      !!----------------------------------------------------------------------
286      This line forced a compilation error
287# endif
288
289#else
290   !!----------------------------------------------------------------------
291   !!   Default option                                 NO LIM sea-ice model
292   !!----------------------------------------------------------------------
293# if defined key_coupled
294      !!----------------------------------------------------------------------
295      !!   'key_coupled' :                            Coupled Ocean/Atmosphere
296      !!----------------------------------------------------------------------
297
298   SUBROUTINE oce_sbc( kt )
299      !!---------------------------------------------------------------------
300      !!                 ***  ROUTINE oce_sbc  ***
301      !!                   
302      !! ** Purpose :   Ocean surface boundaries conditions in
303      !!                coupled ocean/atmosphere case without sea-ice
304      !!
305      !! History :
306      !!   1.0  !  00-10  (O. Marti)  Original code
307      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
308      !!----------------------------------------------------------------------
309      !! * Modules used
310      USE cpl_oce       ! coupled ocean-atmosphere variables
311
312      !! * Arguments
313      INTEGER, INTENT( in  ) ::   kt   ! ocean time step index
314
315      !! * Local declarations
316      INTEGER  ::   ji, jj, jf         ! dummy loop indices
317      REAL(wp) ::   ztgel,          &  ! temporary scalars
318         zice, zhemis, zqrp, zqri,  &  !    "         "
319         zq, zqi, zerp, ze, zei, zro   !    "         "
320      !!----------------------------------------------------------------------
321
322      ! Compute fluxes
323      ! --------------
324
325      DO jj = 1, jpj
326         DO ji = 1, jpi
327
328            ztgel = fzptn(ji,jj)   ! local freezing temperature
329
330            ! opa model ice freeze()
331
332            zice = tmask(ji,jj,1)
333            IF( tn(ji,jj,1) >=  ztgel ) zice = 0.
334            freeze(ji,jj) = zice
335
336            ! hemisphere indicator (=1 north, =-1 south)
337
338            zhemis = float(isign(1, mjg(jj)-(jpjglo/2+1)))
339
340            ! a) net downward radiative flux qsr()
341            ! - AGCM qsrc if no ice
342            ! - zero under ice (zice=1)
343
344            qsr(ji,jj) = (1.-zice)*qsrc(ji,jj)*tmask(ji,jj,1)
345
346            ! b) heat flux damping term qrp()
347            ! - no damping if no  ice      (zice=0)
348            ! - gamma*min(0,t-tgel) if ice (zice=1)
349
350            zqrp = 0.
351            zqri = dqdt0*MIN( 0., tb(ji,jj,1)-ztgel )
352            qrp(ji,jj) = ( ( 1. - zice ) * zqrp + zice * zqri ) * tmask(ji,jj,1)
353
354
355            ! c) net downward heat flux q() = q0 + qrp()
356            ! for q0
357            ! - AGCM qc if no  ice (zice=0)
358            ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if ice (zice=1)
359            zq  = qc(ji,jj)
360            zqi = -3. + zhemis
361            qt(ji,jj) = ( (1.-zice) * zq + zice * zqi ) * tmask(ji,jj,1) + qrp(ji,jj)
362
363            ! d) water flux damping term erp()
364            ! - no damping
365            zerp = 0.
366            erp(ji,jj) = zerp
367
368            ! e) net upward water flux e() = eo + runoff() + erp()
369            ! for e0
370            ! - AGCM if no ice (zice=0)
371            ! - 1.mm/day if climatological and opa ice (zice=1)
372            ze  = ec(ji,jj)
373            zei = 1./rday
374            zro = runoff(ji,jj)
375            emp(ji,jj) = ( ( 1. - zice ) *  ze + zice * zei + zro ) * tmask(ji,jj,1) + erp(ji,jj)
376
377            ! f) net upward water flux for the salinity surface
378            !    boundary condition
379            emps(:,:) = emp(:,:)
380
381         END DO
382      END DO
383
384   END SUBROUTINE oce_sbc
385# elif defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_forced_daily || defined key_flx_core
386      !!-------------------------------------------------------------------------
387      !!   'key_flx_bulk_monthly' or 'key_flx_bulk_daily' or  core   bulk formulea
388      !!   'key_flx_forced_daily'                                or no bulk case
389      !!-------------------------------------------------------------------------
390
391   SUBROUTINE oce_sbc( kt )
392      !!---------------------------------------------------------------------
393      !!                   ***  ROUTINE oce_sbc  ***
394      !!                   
395      !! ** Purpose :   Ocean surface boundary conditions in forced mode
396      !!      using either flux or bulk formulation.
397      !!
398      !! History :
399      !!   1.0  !  99-11  (M. Imbard)  Original code
400      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
401      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
402      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
403      !!----------------------------------------------------------------------
404      !! * Modules used
405      USE daymod                 ! calendar
406#if ! defined key_dtasst
407      USE dtasst, ONLY : rclice  ! sea surface temperature data
408#endif
409#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core
410      USE blk_oce                ! bulk variables
411#endif
412#if defined key_flx_forced_daily
413      USE flx_oce                ! sea-ice/ocean forcings variables
414#endif
415
416      !! * arguments
417      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
418
419      !! * local declarations
420      INTEGER ::   ji, jj        ! dummy loop arguments
421      INTEGER ::   i15, ifreq             !     
422      REAL(wp) ::  zxy
423      REAL(wp) ::  zsice, zqri, zqrp, ztdta, zqrj
424      REAL(wp) ::  zq, zqi, zhemis
425      REAL(wp), DIMENSION(jpi,jpj) :: zeri, zerps, ziclim
426      REAL(wp), DIMENSION(jpi,jpj) :: zqt, zqsr, zemp 
427      !!----------------------------------------------------------------------
428 
429
430#if defined key_flx_core
431      CALL ctl_stop( 'flxcore and no ice model not tested yet' )
432#endif
433
434      ! 1. initialization to zero at kt = nit000
435      ! ---------------------------------------
436
437      IF( kt == nit000 ) THEN     
438         qsr    (:,:) = 0.e0
439         freeze (:,:) = 0.e0
440         qt     (:,:) = 0.e0
441         qla    (:,:) = 0.e0
442         qlw    (:,:) = 0.e0
443         qsb    (:,:) = 0.e0
444         qrp    (:,:) = 0.e0
445         emp    (:,:) = 0.e0
446         emps   (:,:) = 0.e0
447         erp    (:,:) = 0.e0
448#if ! defined key_dynspg_rl
449         dmp    (:,:) = 0.e0
450#endif
451      ENDIF
452
453#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core
454      ifreq      = nfbulk
455      zqt (:,:)  = qsr_oce(:,:) + qnsr_oce(:,:)
456      zqsr(:,:)  = qsr_oce(:,:)
457      zemp(:,:)  = evap(:,:) - tprecip(:,:)
458#endif 
459
460#if defined key_flx_forced_daily
461      ifreq      = 1
462      zqt (:,:)  = p_qt (:,:)
463      zqsr(:,:)  = p_qsr(:,:)
464      zemp(:,:)  = p_emp(:,:)
465#endif 
466#if defined key_flx_core
467      qla(:,:) = (1 - freeze(:,:))*qla_oce + freeze(:,:)*qla_ice(:,:)
468      qsb(:,:) = (1 - freeze(:,:))*qsb_oce + freeze(:,:)*qsb_ice(:,:)
469      qlw(:,:) = (1 - freeze(:,:))*qlw_oce + freeze(:,:)*qlw_ice(:,:)
470#endif
471      IF( MOD( kt-1, ifreq) == 0 ) THEN
472         ! Computation of internal and evaporation damping terms       
473         CALL oce_sbc_dmp
474
475         zsice = - 0.04 / 0.8    ! ratio of isohaline compressibility over isotherme compressibility
476                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
477         ! Flux computation
478         DO jj = 1, jpj
479            DO ji = 1, jpi     
480               ! climatological ice
481               ziclim(ji,jj) = FLOAT( NINT( rclice(ji,jj,1) ) )
482
483               ! avoid surfreezing point           
484               tn(ji,jj,1) = MAX( tn(ji,jj,1), fzptn(ji,jj) )
485
486               ! hemisphere indicator (=1 north, =-1 south)           
487               zhemis = FLOAT( isign(1, mjg(jj) - (jpjdta/2+1) ) )
488
489               ! restoring temperature (ztdta >= to local freezing temperature)           
490#if defined key_dtasst
491               ztdta = MAX( sst(ji,jj),    fzptn(ji,jj) )
492#else
493               ztdta = MAX( t_dta(ji,jj,1), fzptn(ji,jj) )
494#endif
495
496               ! a) net downward radiative flux qsr()           
497               qsr(ji,jj) = (1.-ziclim(ji,jj)) * zqsr(ji,jj) * tmask(ji,jj,1)
498
499               ! b) heat flux damping term qrp()
500               ! - gamma*(t-tlevitus) if no  climatological ice (ziclim=0)
501               ! - gamma*(t-(tgel-1.))  if climatological ice and no opa ice   (ziclim=1 zicopa=0)
502               ! - gamma*min(0,t-tgel) if climatological and opa ice (ziclim=1 zicopa=1)
503
504               zqri = dqdt0 * ( tb(ji,jj,1) - ( fzptn(ji,jj) - 1.) )
505               zqrj = dqdt0 * MIN( 0., tb(ji,jj,1) - fzptn(ji,jj) )
506
507               qrp(ji,jj) =  ( ziclim(ji,jj) * ( (1 - freeze(ji,jj)) * zqri    &
508                 &                                  + freeze(ji,jj)  * zqrj ) ) * tmask(ji,jj,1)
509
510#if ! defined key_flx_bulk_monthly && ! defined key_flx_bulk_daily && ! defined key_flx_core
511               zqrp = dqdt0 * ( tb(ji,jj,1) - ztdta )
512               qrp(ji,jj) = qrp(ji,jj) + (1. - ziclim(ji,jj)) * zqrp
513# endif
514
515               ! c) net downward heat flux q() = q0 + qrp()
516               ! for q0
517               ! - ECMWF fluxes if no climatological ice      (ziclim=0)
518               ! - qrp if climatological ice and no opa ice   (ziclim=1 zicopa=0)
519               ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if climatological and opa ice
520               !                                              (ziclim=1 zicopa=1)
521               zq  = zqt(ji,jj)
522               zqi = -3. + zhemis
523               qt (ji,jj) = ( (1.-ziclim(ji,jj)) * zq   &
524                  +ziclim(ji,jj)  * freeze(ji,jj) * zqi )   &
525                  * tmask(ji,jj,1)   &
526                  + qrp(ji,jj)
527
528            END DO
529         END DO
530
531#if ! defined key_dynspg_rl
532         ! Free-surface
533
534         ! Water flux for zero buoyancy flux if no opa ice and ice clim
535         zeri(:,:) = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
536         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
537
538         ! Contribution to sea level:
539         ! net upward water flux emp() = e-p + runoff() + erp() + dmp + empold
540         emp (:,:) = zemp(:,:)     &   ! e-p data
541            &      + runoff(:,:)   &   ! runoff data
542            &      + erp(:,:)      &   ! restoring term to SSS data
543            &      + dmp(:,:)      &   ! freshwater flux associated with internal damping
544            &      + empold            ! domain averaged annual mean correction
545
546         ! Contribution to salinity:
547         ! net upward water flux emps() = e-p + runoff() + erp() + zerps + empold
548         emps(:,:) = zemp(:,:)     &
549            &      + runoff(:,:)   &
550            &      + erp(:,:)      &
551            &      + zerps(:,:)    &
552            &      + empold 
553#else
554         ! Rigid-lid (emp=emps=E-P-R+Erp)
555         ! freshwater flux
556         zeri(:,:)  = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
557         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
558         emps (:,:) = zemp(:,:)     &
559            &       + runoff(:,:)   &
560            &       + erp(:,:)      &
561            &       + zerps(:,:)
562         emp (:,:) = emps(:,:)
563#endif 
564
565         ! Boundary condition on emp for free surface option
566         ! -------------------------------------------------
567         CALL lbc_lnk( emp, 'T', 1. )
568
569      ENDIF
570
571   END SUBROUTINE oce_sbc
572
573# else
574      !!----------------------------------------------------------------------
575      !!   Default option :                                 Analytical forcing
576      !!----------------------------------------------------------------------
577
578   SUBROUTINE oce_sbc( kt )
579      !!---------------------------------------------------------------------
580      !!                    ***  ROUTINE oce_sbc  ***
581      !!             
582      !! ** Purpose :   provide the thermohaline fluxes (heat and freshwater)
583      !!                to the ocean at each time step.
584      !!
585      !! ** Method  :   Constant surface fluxes (read in namelist (namflx))
586      !!
587      !! ** Action  : - qt, qsr, emp, emps, qrp, erp
588      !!
589      !! History :
590      !!        !  91-03  ()  Original code
591      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
592      !!   9.0  !  04-05  (A. Koch-Larrouy) Add Gyre configuration
593      !!----------------------------------------------------------------------
594      !! * Modules used
595      USE flxrnf                       ! ocean runoffs
596      USE daymod, ONLY : nyear         ! calendar
597      USE dtasss                       ! sea surface salinity data
598
599      !! * arguments
600      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
601
602      !! * local declarations
603      REAL(wp) ::                   & !!! surface fluxes namelist (namflx)
604         q0   = 0.e0,               &  ! net heat flux
605         qsr0 = 0.e0,               &  ! solar heat flux
606         emp0 = 0.e0                   ! net freshwater flux
607      REAL(wp) ::   zsrp,           &
608         zemp_S, zemp_N, zemp_sais, &
609         zTstar, zcos_sais1, zconv, &
610         zcos_sais2
611      REAL(wp) ::           &
612         zsumemp,           &          ! tampon used for the emp sum
613         zsurf,             &          ! tampon used for the domain sum
614         ztime,             &          ! time in hour
615         ztimemax1, ztimemin1, &       ! 21th june,   and 21th december if date0 = 1st january
616         ztimemax2, ztimemin2          ! 21th august, and 21th february if date0 = 1st january
617      REAL(wp), DIMENSION(jpi,jpj) :: t_star
618      INTEGER  ::   ji, jj             ! dummy loop indices
619
620      INTEGER  ::           &
621         zyear0,            &          ! initial year
622         zmonth0,           &          ! initial month
623         zday0,             &          ! initial day
624         zday_year0                    ! initial day since january 1st
625
626      NAMELIST/namflx/ q0, qsr0, emp0
627      !!---------------------------------------------------------------------
628
629      !same temperature, E-P as in HAZELEGER 2000
630
631      IF( cp_cfg == 'gyre' ) THEN
632
633          zyear0  = ndate0 / 10000                              ! initial year
634          zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100           ! initial month
635          zday0   = ndate0 - zyear0 * 10000 - zmonth0 * 100     ! initial day betwen 1 and 30
636          zday_year0= (zmonth0-1)*30.+zday0                     ! initial day betwen 1 and 360
637
638          ! current day (in hours) since january the 1st of the current year
639          ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &    !  total incrementation (in hours)
640             &     - (nyear  - 1) * rjjhh * raajj           !  minus years since beginning of experiment (in hours)
641
642          ! 21th june at 24h  in hours
643          ztimemax1 = ((5.*30.)+21.)* 24.
644          ! 21th december in hours
645          ! rjjhh * raajj / 4 = 1 seasonal cycle in hours
646          ztimemin1 = ztimemax1 + rjjhh * raajj / 2         
647          ! 21th july at 24h in hours
648          ztimemax2 = ((6.*30.)+21.)* 24.
649          ! 21th january day in hours
650          ! rjjhh * raajj / 4 = 1 seasonal cycle in hours
651          ztimemin2 = ztimemax2 - rjjhh * raajj / 2         
652
653          ! amplitudes
654          zemp_S   = 0.7       ! intensity of COS in the South
655          zemp_N   = 0.8       ! intensity of COS in the North
656          zemp_sais= 0.1
657          zTstar   = 28.3      ! intemsity from 28.3 a -5 deg
658
659          ! 1/2 period between 21th June and 21th December
660          zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi )   
661
662          ! 1/2 period between 21th July and 21th January
663          zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
664
665          zconv = 3.16e-5      ! convert 1m/yr->3.16e-5mm/s
666          DO jj = 1, jpj
667             DO ji = 1, jpi
668                ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
669                ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
670                ! 64.5 in summer, 42.5 in winter
671                t_star (ji,jj) = zTstar * ( 1 + 1. / 50. * zcos_sais2 )   &                 
672                   &             * COS( rpi * (gphit(ji,jj) - 5.)        &                 
673                   &                        / (53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) ) 
674                qt  (ji,jj) = dqdt0 * ( tb(ji,jj,1) - t_star(ji,jj) )
675                IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj)) THEN
676                   ! zero at 37.8 deg, max at 24.6 deg
677                   emp  (ji,jj) =   zemp_S * zconv   &
678                      &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  & 
679                      &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
680                   emps (ji,jj) = emp  (ji,jj)
681                ELSE
682                   ! zero at 37.8 deg, max at 46.8 deg
683                   emp (ji,jj) =  - zemp_N * zconv   &
684                      &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  & 
685                      &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
686                   emps (ji,jj) = emp  (ji,jj)
687                ENDIF
688                ! 23.5 deg : tropics
689                qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )   
690             END DO
691          END DO
692
693          ! Compute the emp flux such as its integration on the whole domain at each time is zero
694          IF( nbench /= 1 .AND. nbit_cmp /= 1 ) THEN
695             zsumemp = 0.e0   ;   zsurf = 0.e0
696             DO jj = 1, jpj
697                DO ji = 1, jpi
698                   zsumemp = zsumemp + emp(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
699                   zsurf =  zsurf + tmask(ji,jj,1) * tmask_i(ji,jj)
700                END DO
701             END DO
702
703             IF( lk_mpp )   CALL mpp_sum( zsumemp  )       ! sum over the global domain
704             IF( lk_mpp )   CALL mpp_sum( zsurf    )       ! sum over the global domain
705
706             ! Default GYRE configuration
707             zsumemp = zsumemp / zsurf
708          ELSE
709             ! Benchmark GYRE configuration (to allow the bit to bit comparison between Mpp/Mono case)
710             zsumemp = 0.e0   ;    zsurf = 0.e0
711          ENDIF
712
713          IF( lk_dtasss ) THEN           ! Sea surface salinity damping
714             ! deds0 from (mm/day) => (mm/s)
715             zsrp = deds0 / rday 
716             erp(:,:) = zsrp * ( sss(:,:) - sb(:,:,1) ) / ( sn(:,:,1) + 1.e-20 )
717          ELSE                           ! NO Sea surface salinity damping
718             zsrp = 0.e0
719             erp(:,:) = 0.e0
720          ENDIF
721
722          ! freshwater fluxes terms
723          DO jj = 1, jpj
724             DO ji = 1, jpi
725                emp(ji, jj) =  emp(ji,jj) - zsumemp * tmask(ji,jj,1)
726                emps(ji, jj)=  emp(ji,jj) + erp(ji,jj)
727             END DO
728          END DO
729
730         IF( kt == nit000 .AND. lwp ) THEN
731            WRITE(numout,*)' ocesbc  : analytical formulation for gyre'
732            WRITE(numout,*)' ~~~~~~~ '
733            WRITE(numout,*)'           nyear      = ', nyear
734            WRITE(numout,*)'           ztime      = ',ztime
735            WRITE(numout,*)'           ztimemax1  = ',ztimemax1
736            WRITE(numout,*)'           ztimemin1  = ',ztimemin1
737            WRITE(numout,*)'           ztimemax2  = ',ztimemax2
738            WRITE(numout,*)'           ztimemin2  = ',ztimemin2
739            WRITE(numout,*)'           zyear0     = ', zyear0
740            WRITE(numout,*)'           zmonth0    = ', zmonth0
741            WRITE(numout,*)'           zday0      = ', zday0
742            WRITE(numout,*)'           zday_year0 = ',zday_year0
743            WRITE(numout,*)'           raajj      = ', raajj
744            WRITE(numout,*)'           zemp_S     = ',zemp_S 
745            WRITE(numout,*)'           zemp_N     = ',zemp_N 
746            WRITE(numout,*)'           zemp_sais  = ',zemp_sais
747            WRITE(numout,*)'           zTstar     = ',zTstar 
748            WRITE(numout,*)'           zsumemp    = ',zsumemp
749            WRITE(numout,*)'           zsurf      = ',zsurf 
750            WRITE(numout,*)'           dqdt0      = ',dqdt0
751            WRITE(numout,*)'           deds0      = ',deds0
752            WRITE(numout,*)'           zconv      = ',zconv 
753         ENDIF
754
755      ELSE
756
757         ! Constant surface fluxes
758
759         IF( kt == nit000 ) THEN
760            IF(lwp) THEN
761               WRITE(numout,*)' '
762               WRITE(numout,*)' ocesbc  : Constant surface fluxes read in namelist'
763               WRITE(numout,*)' ~~~~~~~ '
764               WRITE(numout,*)'           Namelist namflx: set the constant flux values'
765               WRITE(numout,*)'              net heat flux          q0   = ', q0  , ' W/m2'
766               WRITE(numout,*)'              solar heat flux        qsr0 = ', qsr0, ' W/m2'
767               WRITE(numout,*)'              net heat flux          emp0 = ', emp0, ' W/m2'
768            ENDIF
769
770            qt    (:,:) = q0
771            qsr   (:,:) = qsr0
772            emp   (:,:) = emp0
773            emps  (:,:) = emp0
774            qrp   (:,:) = 0.e0
775            erp   (:,:) = 0.e0
776
777            runoff(:,:) = 0.e0
778         ENDIF
779
780      ENDIF
781
782   END SUBROUTINE oce_sbc
783
784# endif
785#endif
786
787#if defined key_dtasal
788   !!----------------------------------------------------------------------
789   !!   'key_dtasal'                                          salinity data
790   !!----------------------------------------------------------------------
791   SUBROUTINE oce_sbc_dmp
792      !!---------------------------------------------------------------------
793      !!                   ***  ROUTINE oce_sbc_dmp  ***
794      !!                   
795      !! ** Purpose : Computation of internal and evaporation damping terms
796      !!        for ocean surface boundary conditions
797      !!
798      !! History :
799      !!   9.0  !  04-01  (G. Madec, C. Ethe)  Original code
800      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
801      !!----------------------------------------------------------------------
802      !! * Local declarations
803      INTEGER ::   ji, jj                   ! dummy loop indices
804      REAL(wp), DIMENSION(jpi,jpj)  :: zsss, zfreeze
805      REAL(wp) ::   zerp, zsrp
806      CHARACTER (len=71) :: charout
807#if ! defined key_dynspg_rl
808      REAL(wp) ::   zwei
809      REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj)
810      REAL(wp) ::   zplus, zminus, zadefi
811# if defined key_tradmp
812      INTEGER jk
813      REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp
814# endif
815#endif
816      !!----------------------------------------------------------------------
817
818#if defined key_ice_lim
819      ! sea ice indicator (1 or 0)
820      DO jj = 1, jpj
821         DO ji = 1, jpi
822            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall) )
823         END DO
824      END DO
825      zsss   (:,:) = sss_io(:,:)
826      zfreeze(:,:) = freezn(:,:)
827#else
828      zsss   (:,:) = sb    (:,:,1)
829      zfreeze(:,:) = freeze(:,:)
830#endif
831
832      ! Initialisation
833      ! --------------
834      ! Restoring coefficients on SST and SSS   
835      zsrp = dqdt0 * ro0cpr * rauw   ! (Kg/m2/s)
836
837#if ! defined key_dynspg_rl
838      ! Free-surface
839         
840      ! Internal damping
841# if defined key_tradmp
842      ! Vertical mean of dampind trend (computed in tradmp module)
843      zstrdmp(:,:) = 0.e0
844      DO jk = 1, jpk
845         zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk)
846      END DO
847      ! volume flux associated to internal damping to climatology
848      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 )
849# else
850      dmp(:,:) = 0.e0            ! No internal damping
851# endif
852     
853      !   evaporation damping term ( Surface restoring )
854      zerpplus (:,:) = 0.e0
855      zerpminus(:,:) = 0.e0
856      zplus  =  15. / rday
857      zminus = -15. / rday
858     
859      DO jj = 1, jpj
860         DO ji = 1, jpi
861            zerp = ( 1. - 2.*upsrnfh(ji,jj) ) * zsrp   &
862               & * ( zsss(ji,jj) - s_dta(ji,jj,1) )     &
863               & / ( zsss(ji,jj) + 1.e-20        )
864           
865            zerp = MIN( zerp, zplus  )
866            zerp = MAX( zerp, zminus )
867            erp(ji,jj) = zerp
868            zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 )
869            zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 )
870         END DO
871      END DO
872
873      aplus  = 0.e0
874      aminus = 0.e0
875
876      IF( nbit_cmp == 1) THEN
877
878         IF(ln_ctl)   THEN
879            WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus
880            CALL prt_ctl_info(charout)
881         ENDIF
882         erp(:,:) = 0.e0
883
884      ELSE
885
886         DO jj = 1, jpj
887            DO ji = 1, jpi
888               zwei   = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
889               aplus  = aplus  + zerpplus (ji,jj) * zwei
890               aminus = aminus - zerpminus(ji,jj) * zwei
891            END DO
892         END DO
893         IF( lk_mpp )   CALL mpp_sum( aplus  )   ! sums over the global domain
894         IF( lk_mpp )   CALL mpp_sum( aminus )
895
896         IF(ln_ctl)   THEN
897            WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus
898            CALL prt_ctl_info(charout)
899         ENDIF
900
901         zadefi = MIN( aplus, aminus )
902         IF( zadefi == 0.e0 ) THEN
903            erp(:,:) = 0.e0
904         ELSE
905            erp(:,:) = zadefi * ( zerpplus(:,:) / aplus + zerpminus(:,:) / aminus )
906         ENDIF
907
908      END IF
909#else
910      ! Rigid-lid (emp=emps=E-P-R+Erp)
911
912      erp(:,:) = ( 1. - zfreeze(:,:) ) * zsrp    &   ! surface restoring term
913         &     * ( zsss(:,:) - s_dta(:,:,1) )     &
914         &     / ( zsss(:,:) + 1.e-20      )
915#endif
916
917   END SUBROUTINE oce_sbc_dmp
918
919#else
920   !!----------------------------------------------------------------------
921   !!   Dummy routine                                      NO salinity data
922   !!----------------------------------------------------------------------
923   SUBROUTINE oce_sbc_dmp         ! Dummy routine
924      WRITE(*,*) 'oce_sbc_dmp: you should not have seen that print! error?'
925   END SUBROUTINE oce_sbc_dmp
926#endif
927
928   !!======================================================================
929END MODULE ocesbc
Note: See TracBrowser for help on using the repository browser.