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

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

CT : BUGFIX003 : Compilation error is solved

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