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

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

CT : BUGFIX111 : remove all special characters and lines too long for the GYRE configuration part

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