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

Last change on this file since 417 was 359, checked in by opalod, 18 years ago

nemo_v1_update_033 : RB + CT : Add new surface pressure gradient algorithms

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