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

Last change on this file since 313 was 296, checked in by opalod, 19 years ago

nemo_v1_update_009 : CT : remove all 2D array q(:,:) since they are no more used

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