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

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.6 KB
Line 
1MODULE ocesbc
2   !!======================================================================
3   !!                     ***  MODULE  ocesbc  ***
4   !!                     Ocean surface boundary conditions
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!    oce_sbc : initialization and namelist read
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE oce            ! dynamics and tracers variables
12   USE dom_oce        ! ocean space domain variables
13   USE cpl_oce        ! coupled ocean-atmosphere variables
14   USE ice_oce
15   USE blk_oce
16   USE phycst         ! Define parameters for the routines
17   USE taumod
18   USE flxmod
19   USE flxrnf
20   USE tradmp         ! damping salinity trend
21   USE dtatem
22   USE dtasal
23   USE ocfzpt
24   USE lbclnk
25   USE in_out_manager ! I/O manager
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Accessibility
31   PUBLIC oce_sbc    ! routine called by step
32
33   !! * Shared module variables
34   REAL(wp), PUBLIC ::   &
35      aplus, aminus,     &
36      empold = 0.e0  ! current year freshwater budget correction
37   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
38      qt  ,         &  ! total surface heat flux (w/m2)
39      q   ,         &  ! surface heat flux (w/m2)
40      qsr ,         &  ! solar radiation (w/m2)
41      emp ,         &  ! evaporation minus precipitation (kg/m2/s = mm/s)
42      emps,         &  ! evaporation - precipitation (free surface)
43      qrp ,         &  ! heat flux damping (w/m2)
44      erp              ! evaporation damping (kg/m2/s = mm/s)
45#if defined key_dynspg_fsc
46   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
47      dmp              ! internal dampind term
48#endif
49
50#  include "domzgr_substitute.h90"
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58#if defined key_ice_lim
59   !!----------------------------------------------------------------------
60   !!   'key_ice_lim' :                                   LIM sea-ice model
61   !!----------------------------------------------------------------------
62# if defined key_coupled
63      !!----------------------------------------------------------------------
64      !!   'key_coupled' :                            Coupled Ocean/Atmosphere
65      !!----------------------------------------------------------------------
66
67   SUBROUTINE oce_sbc( kt )
68      !!---------------------------------------------------------------------
69      !!                 ***  ROUTINE oce_sbc  ***
70      !!                   
71      !! ** Purpose :   Ocean surface boundaries conditions with
72      !!        Louvain la Neuve Sea Ice Model in coupled mode
73      !!
74      !! History :
75      !!   1.0  !  00-10  (O. Marti)  Original code
76      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
77      !!----------------------------------------------------------------------
78      !! * Local declarations
79      INTEGER ::   ji, jj                   ! dummy loop indices
80      REAL(wp) ::   ztx, ztaux, zty, ztauy, ztrp, zsrp
81      REAL(wp) ::   ztdta, ztgel, zqrp
82      !!----------------------------------------------------------------------
83 
84      ! 1. initialization to zero at kt = nit000
85      ! ---------------------------------------
86     
87      IF( kt == nit000 ) THEN     
88         qsr   (:,:) = 0.e0
89         freeze(:,:) = 0.e0
90         qt    (:,:) = 0.e0
91         q     (:,:) = 0.e0
92         qrp   (:,:) = 0.e0
93         emp   (:,:) = 0.e0
94         emps  (:,:) = 0.e0
95         erp   (:,:) = 0.e0
96#if defined key_dynspg_fsc 
97         dmp   (:,:) = 0.e0
98#endif
99      ENDIF
100
101      IF( MOD( kt-1, nfice ) == 0 ) THEN
102
103         CALL oce_sbc_dmp   ! Computation of internal and evaporation damping terms       
104
105         ! Surface heat flux (W/m2)
106         ! -----------------------
107
108         ! restoring heat flux
109         DO jj = 1, jpj
110            DO ji = 1, jpi
111               ztgel = fzptn(ji,jj)
112#if defined key_dtasst
113               ztdta = MAX( sst(ji,jj),    ztgel )
114#else
115               ztdta = MAX( t_dta(ji,jj,1), ztgel )
116#endif
117               zqrp = ztrp * ( tb(ji,jj,1) - ztdta )
118
119               qrp(ji,jj) = (1.0-freeze(ji,jj) ) * zqrp
120            END DO
121         END DO
122
123         ! non solar heat flux + solar flux + restoring
124         q  (:,:) = fnsolar(:,:) + fsolar(:,:) + qrp(:,:)
125         qt (:,:) = q(:,:)
126
127         ! solar flux
128         qsr(:,:) = fsolar(:,:)
129
130#if defined key_dynspg_fsc 
131         ! total concentration/dilution effect (use on SSS)
132         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
133
134         ! total volume flux (use on sea-surface height)
135         emp (:,:) = fmass(:,:)  -  dmp(:,:) + runoff(:,:) + erp(:,:)
136#else
137         ! Rigid-lid (emp=emps=E-P-R+Erp)
138         ! freshwater flux
139         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
140         emp (:,:) = emps(:,:)
141#endif
142
143         DO jj = 1, jpjm1
144            DO ji = 1, fs_jpim1   ! vertor opt.
145               ztx   = 0.5 * ( freeze(ji+1,jj) + freeze(ji+1,jj+1) )
146               ztaux = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) )
147               taux(ji,jj) = (1.0-ztx) * taux(ji,jj) + ztx * ztaux
148
149               zty   = 0.5 * ( freeze(ji,jj+1) + freeze(ji+1,jj+1) )
150               ztauy = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) )
151               tauy(ji,jj) = (1.0-zty) * tauy(ji,jj) + zty * ztauy
152            END DO
153         END DO
154         CALL lbc_lnk( taux, 'U', -1. )
155         CALL lbc_lnk( tauy, 'V', -1. )   
156
157         ! Re-initialization of fluxes
158         sst_io(:,:) = 0.0
159         sss_io(:,:) = 0.0
160         u_io  (:,:) = 0.0
161         v_io  (:,:) = 0.0
162         gtaux (:,:) = 0.
163         gtauy (:,:) = 0.
164
165      ENDIF
166
167   END SUBROUTINE oce_sbc
168
169# elif defined key_flx_bulk_monthly || defined key_flx_bulk_daily
170      !!----------------------------------------------------------------------
171      !!   'key_ice_lim'                              with  LIM sea-ice model
172      !!----------------------------------------------------------------------
173
174   SUBROUTINE oce_sbc( kt )
175      !!---------------------------------------------------------------------
176      !!                   ***  ROUTINE oce_sbc  ***
177      !!                   
178      !! ** Purpose : - Ocean surface boundary conditions with LIM sea-ice
179      !!        model in forced mode using bulk formulea
180      !!
181      !! History :
182      !!   1.0  !  99-11  (M. Imbard)  Original code
183      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
184      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
185      !!----------------------------------------------------------------------
186      !! * arguments
187      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
188
189      !! * Local declarations
190      INTEGER ::   ji, jj                   ! dummy loop indices
191      REAL ztx, ztaux, zty, ztauy
192      !!----------------------------------------------------------------------
193
194      ! 1. initialization to zero at kt = nit000
195      ! ---------------------------------------
196     
197      IF( kt == nit000 ) THEN     
198         qsr    (:,:) = 0.e0
199         qt     (:,:) = 0.e0
200         q      (:,:) = 0.e0
201         qrp    (:,:) = 0.e0
202         emp    (:,:) = 0.e0
203         emps   (:,:) = 0.e0
204         erp    (:,:) = 0.e0
205#if defined key_dynspg_fsc 
206         dmp    (:,:) = 0.e0
207#endif
208      ENDIF
209
210      IF( MOD( kt-1, nfice ) == 0 ) THEN
211
212         CALL oce_sbc_dmp       ! Computation of internal and evaporation damping terms       
213
214         ! Surface Ocean fluxes
215         ! ====================
216
217         ! Surface heat flux (W/m2)
218         ! -----------------
219
220         q   (:,:) = fnsolar(:,:) + fsolar(:,:)     ! non solar heat flux + solar flux
221         qt  (:,:) = q(:,:)
222         qsr (:,:) = fsolar(:,:)                     ! solar flux
223
224#if defined key_dynspg_fsc     
225         ! total concentration/dilution effect (use on SSS)
226         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) + empold
227
228         ! total volume flux (use on sea-surface height)
229         emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold     
230#else
231         ! Rigid-lid (emp=emps=E-P-R+Erp)
232         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)     ! freshwater flux
233         emp (:,:) = emps(:,:)
234
235#endif
236
237         ! Surface stress
238         ! --------------
239
240         ! update the stress beloww sea-ice area
241         DO jj = 1, jpjm1
242            DO ji = 1, fs_jpim1   ! vertor opt.
243               ztx         = MAX( freezn(ji,jj), freezn(ji,jj+1) )   ! ice/ocean indicator at U- and V-points
244               zty         = MAX( freezn(ji,jj), freezn(ji+1,jj) )
245               ztaux       = 0.5 *( ftaux(ji+1,jj) + ftaux(ji+1,jj+1) ) ! ice-ocean stress at U- and V-points
246               ztauy       = 0.5 *( ftauy(ji,jj+1) + ftauy(ji+1,jj+1) )
247               taux(ji,jj) = (1.-ztx) * taux(ji,jj) + ztx * ztaux    ! stress at the ocean surface
248               tauy(ji,jj) = (1.-zty) * tauy(ji,jj) + zty * ztauy
249            END DO
250         END DO
251
252         ! boundary condition on the stress (taux,tauy)
253         CALL lbc_lnk( taux, 'U', -1. )
254         CALL lbc_lnk( tauy, 'V', -1. )
255
256         ! Re-initialization of fluxes
257         sst_io(:,:) = 0.0
258         sss_io(:,:) = 0.0
259         u_io  (:,:) = 0.0
260         v_io  (:,:) = 0.0
261
262      ENDIF
263
264   END SUBROUTINE oce_sbc
265
266# else
267      !!----------------------------------------------------------------------
268      !!   Error option               LIM sea-ice model requires bulk formulea
269      !!----------------------------------------------------------------------
270      This line forced a compilation error
271# endif
272
273#else
274   !!----------------------------------------------------------------------
275   !!   Default option                                 NO LIM sea-ice model
276   !!----------------------------------------------------------------------
277# if defined key_coupled
278      !!----------------------------------------------------------------------
279      !!   'key_coupled' :                            Coupled Ocean/Atmosphere
280      !!----------------------------------------------------------------------
281
282   SUBROUTINE oce_sbc( kt )
283      !!---------------------------------------------------------------------
284      !!                 ***  ROUTINE oce_sbc  ***
285      !!                   
286      !! ** Purpose :   Ocean surface boundaries conditions in
287      !!                coupled ocean/atmosphere case without sea-ice
288      !!
289      !! History :
290      !!   1.0  !  00-10  (O. Marti)  Original code
291      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
292      !!----------------------------------------------------------------------
293      !! * Modules used
294      USE cpl_oce 
295
296      !! * Arguments
297      INTEGER, INTENT( in  ) ::   kt   ! ocean time step index
298
299      !! * Local declarations
300      INTEGER  ::   ji, jj, jf         ! dummy loop indices
301      REAL(wp) ::   ztrp, ztgel,    &  ! temporary scalars
302         zice, zhemis, zqrp, zqri,  &  !    "         "
303         zq, zqi, zerp, ze, zei, zro   !    "         "
304      !!----------------------------------------------------------------------
305
306      ! Compute fluxes
307      ! --------------
308
309      ! constant initialization
310      ztrp = -40.   ! restoring term for temperature (w/m2/k)
311
312      DO jj = 1, jpj
313         DO ji = 1, jpi
314
315            ztgel = fzptn(ji,jj)   ! local freezing temperature
316
317            ! opa model ice freeze()
318
319            zice = tmask(ji,jj,1)
320            IF( tn(ji,jj,1) >=  ztgel ) zice = 0.
321            freeze(ji,jj) = zice
322
323            ! hemisphere indicator (=1 north, =-1 south)
324           
325            zhemis = float(isign(1, mjg(jj)-(jpjglo/2+1)))
326           
327            ! a) net downward radiative flux qsr()
328            ! - AGCM qsrc if no ice
329            ! - zero under ice (zice=1)
330
331            qsr(ji,jj) = (1.-zice)*qsrc(ji,jj)*tmask(ji,jj,1)
332
333            ! b) heat flux damping term qrp()
334            ! - no damping if no  ice      (zice=0)
335            ! - gamma*min(0,t-tgel) if ice (zice=1)
336
337            zqrp = 0.
338            zqri = ztrp*MIN( 0., tb(ji,jj,1)-ztgel )
339            qrp(ji,jj) = ( ( 1. - zice ) * zqrp + zice * zqri ) * tmask(ji,jj,1)
340
341
342            ! c) net downward heat flux q() = q0 + qrp()
343            ! for q0
344            ! - AGCM qc if no  ice (zice=0)
345            ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if ice (zice=1)
346            zq  = qc(ji,jj)
347            zqi = -3. + zhemis
348            qt(ji,jj) = ( (1.-zice) * zq + zice * zqi ) * tmask(ji,jj,1) + qrp(ji,jj)
349            q (ji,jj) = qt(ji,jj)
350           
351            ! d) water flux damping term erp()
352            ! - no damping
353            zerp = 0.
354            erp(ji,jj) = zerp
355           
356            ! e) net upward water flux e() = eo + runoff() + erp()
357            ! for e0
358            ! - AGCM if no ice (zice=0)
359            ! - 1.mm/day if climatological and opa ice (zice=1)
360            ze  = ec(ji,jj)
361            zei = 1./rday
362            zro = runoff(ji,jj)
363            emp(ji,jj) = ( ( 1. - zice ) *  ze + zice * zei + zro ) * tmask(ji,jj,1) + erp(ji,jj)
364           
365         END DO
366      END DO
367
368   END SUBROUTINE oce_sbc
369
370# elif defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_forced_daily
371      !!-------------------------------------------------------------------------
372      !!   'key_flx_bulk_monthly' or 'key_flx_bulk_daily' or        bulk formulea
373      !!   'key_flx_forced_daily'                                or no bulk case
374      !!-------------------------------------------------------------------------
375
376   SUBROUTINE oce_sbc( kt )
377      !!---------------------------------------------------------------------
378      !!                   ***  ROUTINE oce_sbc  ***
379      !!                   
380      !! ** Purpose :   Ocean surface boundary conditions in forced mode
381      !!      using either flux or bulk formulation.
382      !!
383      !! History :
384      !!   1.0  !  99-11  (M. Imbard)  Original code
385      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
386      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
387      !!----------------------------------------------------------------------
388      !! * Modules used
389      USE daymod
390#if ! defined key_dtasst
391      USE dtasst, ONLY : rclice
392#endif
393#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
394      USE blk_oce
395#endif
396#if defined key_flx_forced_daily
397      USE flx_oce
398#endif
399
400      !! * arguments
401      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
402
403      !! * local declarations
404      INTEGER ::   ji, jj        ! dummy loop arguments
405      INTEGER ::   i15, ifreq             !     
406      REAL(wp) ::  zxy
407      REAL(wp) ::  zsice, zqri, zqrp, ztdta, zqrj
408      REAL(wp) ::  zq, zqi, zhemis, ztrp
409      REAL(wp), DIMENSION(jpi,jpj) :: zeri, zerps, ziclim
410      REAL(wp), DIMENSION(jpi,jpj) :: zqt, zqsr, zemp 
411      !!----------------------------------------------------------------------
412 
413      ! 1. initialization to zero at kt = nit000
414      ! ---------------------------------------
415     
416      IF( kt == nit000 ) THEN     
417         qsr    (:,:) = 0.e0
418         freeze (:,:) = 0.e0
419         qt     (:,:) = 0.e0
420         q      (:,:) = 0.e0
421         qrp    (:,:) = 0.e0
422         emp    (:,:) = 0.e0
423         emps   (:,:) = 0.e0
424         erp    (:,:) = 0.e0
425#if defined key_dynspg_fsc 
426         dmp    (:,:) = 0.e0
427#endif
428      ENDIF
429
430#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
431      ifreq      = nfbulk
432      zqt (:,:)  = qsr_oce(:,:) + qnsr_oce(:,:)
433      zqsr(:,:)  = qsr_oce(:,:)
434      zemp(:,:)  = evap(:,:) - tprecip(:,:)
435#endif 
436
437#if defined key_flx_forced_daily
438      ifreq      = 1
439      zqt (:,:)  = p_qt (:,:)
440      zqsr(:,:)  = p_qsr(:,:)
441      zemp(:,:)  = p_emp(:,:)
442#endif
443
444      IF( MOD( kt-1, ifreq) == 0 ) THEN
445         ! Computation of internal and evaporation damping terms       
446         CALL oce_sbc_dmp
447
448         ztrp = -40.             ! restoring terme for temperature (w/m2/k)   
449         zsice = - 0.04 / 0.8    ! ratio of isohaline compressibility over isotherme compressibility
450                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
451         ! Flux computation
452         DO jj = 1, jpj
453            DO ji = 1, jpi     
454               ! climatological ice
455#if defined key_dtasst 
456               ziclim(ji,jj) = FLOAT( NINT( rclice(ji,jj,1) ) )
457#else
458               ! tested only with key key_dtasst (A. Lazar 07/2001)
459               ! this loop in CASE of interpolation of monthly rclice
460               i15           = INT( 2.* FLOAT(nday) / (FLOAT( nobis(nmonth) ) + 0.5) )
461               zxy           = FLOAT(nday) / FLOAT(nobis(nmonth)) + 0.5 - i15
462               ziclim(ji,jj) = FLOAT( NINT( (1-zxy) * rclice(ji,jj,1) + zxy  * rclice(ji,jj,2) ) )
463#endif
464
465               ! avoid surfreezing point           
466               tn(ji,jj,1) = MAX( tn(ji,jj,1), fzptn(ji,jj) )
467
468               ! hemisphere indicator (=1 north, =-1 south)           
469               zhemis = FLOAT( isign(1, mjg(jj) - (jpjdta/2+1) ) )
470
471               ! restoring temperature (ztdta >= to local freezing temperature)           
472#if defined key_dtasst
473               ztdta = MAX( sst(ji,jj),    fzptn(ji,jj) )
474#else
475               ztdta = MAX( t_dta(ji,jj,1), fzptn(ji,jj) )
476#endif
477
478               ! a) net downward radiative flux qsr()           
479               qsr(ji,jj) = zqsr(ji,jj) * tmask(ji,jj,1)
480
481               ! b) heat flux damping term qrp()
482               ! - gamma*(t-tlevitus) if no  climatological ice (ziclim=0)
483               ! - gamma*(t-(tgel-1.))  if climatological ice and no opa ice   (ziclim=1 zicopa=0)
484               ! - gamma*min(0,t-tgel) if climatological and opa ice (ziclim=1 zicopa=1)
485
486               zqrp = ztrp * ( tb(ji,jj,1) - ztdta )
487               zqri = ztrp * ( tb(ji,jj,1) - ( fzptn(ji,jj) - 1.) )
488               zqrj = ztrp * MIN( 0., tb(ji,jj,1) - fzptn(ji,jj) )
489               qrp(ji,jj) = ( (1. - ziclim(ji,jj)) * zqrp   &
490                  + ziclim(ji,jj)  * ( ( 1 - freeze(ji,jj)) * zqri   &
491                  + freeze(ji,jj)  * zqrj ) ) * tmask(ji,jj,1)
492
493               ! c) net downward heat flux q() = q0 + qrp()
494               ! for q0
495               ! - ECMWF fluxes if no climatological ice      (ziclim=0)
496               ! - qrp if climatological ice and no opa ice   (ziclim=1 zicopa=0)
497               ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if climatological and opa ice
498               !                                              (ziclim=1 zicopa=1)
499               zq  = zqt(ji,jj)
500               zqi = -3. + zhemis
501               qt (ji,jj) = ( (1.-ziclim(ji,jj)) * zq   &
502                  +ziclim(ji,jj)  * freeze(ji,jj) * zqi )   &
503                  * tmask(ji,jj,1)   &
504                  + qrp(ji,jj)
505               q  (ji,jj) = qt (ji,jj)
506
507            END DO
508         END DO
509
510#if defined key_dynspg_fsc 
511         ! Free-surface
512
513         ! Water flux for zero buoyancy flux if no opa ice and ice clim
514         zeri(:,:) = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
515         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
516
517         ! Contribution to sea level:
518         ! net upward water flux emp() = e-p + runoff() + erp() + dmp + empold
519         emp (:,:) = zemp(:,:)     &   ! e-p data
520            &      + runoff(:,:)   &   ! runoff data
521            &      + erp(:,:)      &   ! restoring term to SSS data
522            &      + dmp(:,:)      &   ! freshwater flux associated with internal damping
523            &      + empold            ! domain averaged annual mean correction
524
525         ! Contribution to salinity:
526         ! net upward water flux emps() = e-p + runoff() + erp() + zerps + empold
527         emps(:,:) = zemp(:,:)     &
528            &      + runoff(:,:)   &
529            &      + erp(:,:)      &
530            &      + zerps(:,:)    &
531            &      + empold 
532#else
533         ! Rigid-lid (emp=emps=E-P-R+Erp)
534         ! freshwater flux
535         zeri(:,:)  = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
536         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
537         emps (:,:) = zemp(:,:)     &
538            &       + runoff(:,:)   &
539            &       + erp(:,:)      &
540            &       + zerps(:,:)
541         emp (:,:) = emps(:,:)
542#endif 
543
544         ! Boundary condition on emp for free surface option
545         ! -------------------------------------------------
546         CALL lbc_lnk( emp, 'T', 1. )
547     
548      ENDIF
549
550   END SUBROUTINE oce_sbc
551
552# else
553      !!----------------------------------------------------------------------
554      !!   Default option :                                 Analytical forcing
555      !!----------------------------------------------------------------------
556
557   SUBROUTINE oce_sbc( kt )
558      !!---------------------------------------------------------------------
559      !!                    ***  ROUTINE oce_sbc  ***
560      !!             
561      !! ** Purpose :   provide the thermohaline fluxes (heat and freshwater)
562      !!      to the ocean at each time step.
563      !!
564      !! ** Method  :   Constant surface fluxes (read in namelist (namflx))
565      !!
566      !! ** Action  : - q, qt, qsr, emp, emps, qrp, erp
567      !!
568      !! History :
569      !!        !  91-03  ()  Original code
570      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
571      !!----------------------------------------------------------------------
572      !! * Modules used
573      USE flxrnf                       ! ocean runoffs
574
575      !! * arguments
576      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
577
578      !! * local declarations
579      REAL(wp) ::                   & !!! surface fluxes namelist (namflx)
580         q0   = 0.e0,               &  ! net heat flux
581         qsr0 = 0.e0,               &  ! solar heat flux
582         emp0 = 0.e0                   ! net freshwater flux
583
584      NAMELIST/namflx/ q0, qsr0, emp0
585      !!---------------------------------------------------------------------
586
587      IF( kt == nit000 ) THEN
588
589         ! Read Namelist namflx : surface thermohaline fluxes
590         ! --------------------
591         REWIND ( numnam )
592         READ   ( numnam, namflx )
593
594         IF(lwp) WRITE(numout,*)' '
595         IF(lwp) WRITE(numout,*)' ocesbc  : Constant surface fluxes read in namelist'
596         IF(lwp) WRITE(numout,*)' ~~~~~~~ '
597         IF(lwp) WRITE(numout,*)'           Namelist namflx: set the constant flux values'
598         IF(lwp) WRITE(numout,*)'              net heat flux          q0   = ', q0  , ' W/m2'
599         IF(lwp) WRITE(numout,*)'              solar heat flux        qsr0 = ', qsr0, ' W/m2'
600         IF(lwp) WRITE(numout,*)'              net heat flux          emp0 = ', emp0, ' W/m2'
601
602         qt    (:,:) = q0
603         qsr   (:,:) = qsr0
604         q     (:,:) = q0
605         emp   (:,:) = emp0
606         emps  (:,:) = emp0
607         qrp   (:,:) = 0.e0
608         erp   (:,:) = 0.e0
609
610         runoff(:,:) = 0.e0
611     ENDIF
612
613  END SUBROUTINE oce_sbc
614
615# endif
616#endif
617
618   SUBROUTINE oce_sbc_dmp
619      !!---------------------------------------------------------------------
620      !!                   ***    ROUTINE oce_sbc_dmp ***
621      !!                   
622      !! ** Purpose : Computation of internal and evaporation damping terms
623      !!        for ocean surface boundary conditions
624      !!
625      !! History :
626      !!   1.0  !  99-11  (M. Imbard)  Original code
627      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
628      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
629      !!----------------------------------------------------------------------
630      !! * Local declarations
631      INTEGER ::   ji, jj                   ! dummy loop indices
632      REAL(wp), DIMENSION(jpi,jpj)  :: zsss, zfreeze
633#if defined key_dynspg_fsc
634      REAL zwei, zerp, ztrp, zsrp
635      REAL zerpplus(jpi,jpj), zerpminus(jpi,jpj)
636      REAL zplus, zminus, zadefi
637# if defined key_tradmp
638      INTEGER jk
639      REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp
640# endif
641#endif
642      !!----------------------------------------------------------------------
643
644#if defined key_ice_lim
645      ! sea ice indicator (1 or 0)
646      DO jj = 1, jpj
647         DO ji = 1, jpi
648            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall))
649         END DO
650      END DO
651      zsss   (:,:) = sss_io(:,:)
652      zfreeze(:,:) = freezn(:,:)
653#else
654      zsss   (:,:) = sn    (:,:,1)
655      zfreeze(:,:) = freeze(:,:)
656#endif
657
658      ! Initialisation
659      ! --------------
660      ! Restoring coefficients on SST and SSS   
661      IF( lk_cpl ) THEN
662         ztrp = 0.e0
663         zsrp = 0.e0
664      ELSE
665         ztrp  = -40.                   ! (W/m2/K)   
666         zsrp  = ztrp * ro0cpr * rauw   ! (Kg/m2/s2)
667      ENDIF
668
669#if defined key_dynspg_fsc 
670      ! Free-surface
671         
672      ! Internal damping
673# if defined key_tradmp
674      ! Vertical mean of dampind trend (computed in tradmp module)
675      zstrdmp(:,:) = 0.e0
676      DO jk = 1, jpk
677         zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk)
678      END DO
679      ! volume flux associated to internal damping to climatology
680      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + rsmall )
681# else
682      dmp(:,:) = 0.e0            ! No internal damping
683# endif
684     
685      !   evaporation damping term ( Surface restoring )
686      zerpplus (:,:) = 0.e0
687      zerpminus(:,:) = 0.e0
688      zplus  =  15. / rday
689      zminus = -15. / rday
690     
691      DO jj = 1, jpj
692         DO ji = 1, jpi
693            zerp = ( 1. - 2.*upsrnfh(ji,jj) ) * zsrp   &
694               & * ( zsss(ji,jj) - s_dta(ji,jj,1) )     &
695               & / ( zsss(ji,jj) + rsmall        )
696           
697            zerp = MIN( zerp, zplus  )
698            zerp = MAX( zerp, zminus )
699            erp(ji,jj) = zerp
700            zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 )
701            zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 )
702         END DO
703      END DO
704
705      aplus  = 0.e0
706      aminus = 0.e0
707      DO jj = 1, jpj
708         DO ji = 1, jpi
709            zwei   = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
710            aplus  = aplus  + zerpplus (ji,jj) * zwei
711            aminus = aminus - zerpminus(ji,jj) * zwei
712         END DO
713      END DO
714# if defined key_mpp
715      CALL mpp_sum( aplus  )   ! mpp: sum over all the global domain
716      CALL mpp_sum( aminus )
717# endif
718      IF( l_ctl .AND. lwp ) WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
719
720      zadefi = MIN( aplus, aminus )
721      IF( zadefi == 0.0 ) THEN
722         erp(:,:) = 0.e0
723      ELSE
724         erp(:,:) = zadefi * ( zerpplus(:,:) / aplus + zerpminus(:,:) / aminus )
725      ENDIF
726#else
727      ! Rigid-lid (emp=emps=E-P-R+Erp)
728     
729      erp(:,:) = ( 1. - zfreeze(:,:) ) * zsrp    &   ! surface restoring term
730         &     * ( zsss(:,:) - s_dta(:,:,1) )     &
731         &     / ( zsss(:,:) + rsmall      )
732#endif
733
734   END SUBROUTINE oce_sbc_dmp
735
736   !!======================================================================
737END MODULE ocesbc
Note: See TracBrowser for help on using the repository browser.