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.
trabbc.F90 in branches/UKMO/dev_r5518_GO6_package_for_static_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_package_for_static_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90 @ 7599

Last change on this file since 7599 was 7599, checked in by timgraham, 8 years ago

Fully working version
File size: 9.9 KB
Line 
1MODULE trabbc
2   !!==============================================================================
3   !!                       ***  MODULE  trabbc  ***
4   !! Ocean active tracers:  bottom boundary condition (geothermal heat flux)
5   !!==============================================================================
6   !! History :  OPA  ! 1999-10 (G. Madec)  original code
7   !!   NEMO     1.0  ! 2002-08 (G. Madec)  free form + modules
8   !!             -   ! 2002-11 (A. Bozec)  tra_bbc_init: original code
9   !!            3.3  ! 2010-10 (G. Madec)  dynamical allocation + suppression of key_trabbc
10   !!             -   ! 2010-11 (G. Madec)  use mbkt array (deepest ocean t-level)
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   tra_bbc      : update the tracer trend at ocean bottom
15   !!   tra_bbc_init : initialization of geothermal heat flux trend
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean variables
18   USE dom_oce         ! domain: ocean
19   USE phycst          ! physical constants
20   USE trd_oce         ! trends: ocean variables
21   USE trdtra          ! trends manager: tracers
22   USE in_out_manager  ! I/O manager
23   USE iom             ! I/O manager
24   USE fldread         ! read input fields
25   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
26   USE lib_mpp           ! distributed memory computing library
27   USE prtctl          ! Print control
28   USE wrk_nemo        ! Memory Allocation
29   USE timing          ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC tra_bbc          ! routine called by step.F90
35   PUBLIC tra_bbc_init     ! routine called by opa.F90
36
37   !                                 !!* Namelist nambbc: bottom boundary condition *
38   LOGICAL, PUBLIC ::   ln_trabbc     !: Geothermal heat flux flag
39   INTEGER         ::   nn_geoflx     !  Geothermal flux (=1:constant flux, =2:read in file )
40   REAL(wp)        ::   rn_geoflx_cst !  Constant value of geothermal heat flux
41
42   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd0   ! geothermal heating trend
43   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read)
44 
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE tra_bbc( kt )
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE tra_bbc  ***
57      !!
58      !! ** Purpose :   Compute the bottom boundary contition on temperature
59      !!              associated with geothermal heating and add it to the
60      !!              general trend of temperature equations.
61      !!
62      !! ** Method  :   The geothermal heat flux set to its constant value of
63      !!              86.4 mW/m2 (Stein and Stein 1992, Huang 1999).
64      !!       The temperature trend associated to this heat flux through the
65      !!       ocean bottom can be computed once and is added to the temperature
66      !!       trend juste above the bottom at each time step:
67      !!            ta = ta + Qsf / (rau0 rcp e3T) for k= mbkt
68      !!       Where Qsf is the geothermal heat flux.
69      !!
70      !! ** Action  : - update the temperature trends (ta) with the trend of
71      !!                the ocean bottom boundary condition
72      !!
73      !! References : Stein, C. A., and S. Stein, 1992, Nature, 359, 123-129.
74      !!              Emile-Geay and Madec, 2009, Ocean Science.
75      !!----------------------------------------------------------------------
76      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
77      !!
78      INTEGER  ::   ji, jj, ik    ! dummy loop indices
79      REAL(wp) ::   zqgh_trd      ! geothermal heat flux trend
80      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt
81      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmp
82      !!----------------------------------------------------------------------
83      !
84      IF( nn_timing == 1 )  CALL timing_start('tra_bbc')
85      !
86      IF( l_trdtra )   THEN         ! Save ta and sa trends
87         CALL wrk_alloc( jpi, jpj, jpk, ztrdt )
88         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
89      ENDIF
90      !
91      IF( kt == nit000 ) THEN
92      !Output the geothermal heat flux once for CMIP6 diagnostics
93         CALL wrk_alloc(jpi, jpj, ztmp)
94         DO jj=1,jpj
95            DO ji=1,jpi
96               ztmp(ji,jj) = qgh_trd0(ji,jj)*rau0_rcp
97            ENDDO
98         ENDDO
99         CALL iom_put( 'hfgeou', ztmp )
100         CALL wrk_dealloc(jpi, jpj, ztmp)
101      ENDIF
102
103      !                             !  Add the geothermal heat flux trend on temperature
104      DO jj = 2, jpjm1
105         DO ji = 2, jpim1
106            ik = mbkt(ji,jj)
107            zqgh_trd = qgh_trd0(ji,jj) / fse3t(ji,jj,ik)
108            tsa(ji,jj,ik,jp_tem) = tsa(ji,jj,ik,jp_tem) + zqgh_trd
109         END DO
110      END DO
111      !
112      CALL lbc_lnk( tsa(:,:,:,jp_tem) , 'T', 1. )
113      !
114      IF( l_trdtra ) THEN        ! Save the geothermal heat flux trend for diagnostics
115         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
116         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt )
117         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )
118      ENDIF
119      !
120      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
121      !
122      IF( nn_timing == 1 )  CALL timing_stop('tra_bbc')
123      !
124   END SUBROUTINE tra_bbc
125
126
127   SUBROUTINE tra_bbc_init
128      !!----------------------------------------------------------------------
129      !!                  ***  ROUTINE tra_bbc_init  ***
130      !!
131      !! ** Purpose :   Compute once for all the trend associated with geothermal
132      !!              heating that will be applied at each time step at the
133      !!              last ocean level
134      !!
135      !! ** Method  :   Read the nambbc namelist and check the parameters.
136      !!
137      !! ** Input   : - Namlist nambbc
138      !!              - NetCDF file  : geothermal_heating.nc ( if necessary )
139      !!
140      !! ** Action  : - read/fix the geothermal heat qgh_trd0
141      !!----------------------------------------------------------------------
142      USE iom
143      !!
144      INTEGER  ::   ji, jj              ! dummy loop indices
145      INTEGER  ::   inum                ! temporary logical unit
146      INTEGER  ::   ios                 ! Local integer output status for namelist read
147      INTEGER  ::   ierror              ! local integer
148      !
149      TYPE(FLD_N)        ::   sn_qgh    ! informations about the geotherm. field to be read
150      CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files
151      !
152      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 
153      !!----------------------------------------------------------------------
154
155      REWIND( numnam_ref )              ! Namelist nambbc in reference namelist : Bottom momentum boundary condition
156      READ  ( numnam_ref, nambbc, IOSTAT = ios, ERR = 901)
157901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbc in reference namelist', lwp )
158
159      REWIND( numnam_cfg )              ! Namelist nambbc in configuration namelist : Bottom momentum boundary condition
160      READ  ( numnam_cfg, nambbc, IOSTAT = ios, ERR = 902 )
161902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbc in configuration namelist', lwp )
162      IF(lwm) WRITE ( numond, nambbc )
163
164      IF(lwp) THEN                     ! Control print
165         WRITE(numout,*)
166         WRITE(numout,*) 'tra_bbc : Bottom Boundary Condition (bbc), apply a Geothermal heating'
167         WRITE(numout,*) '~~~~~~~   '
168         WRITE(numout,*) '   Namelist nambbc : set bbc parameters'
169         WRITE(numout,*) '      Apply a geothermal heating at ocean bottom   ln_trabbc     = ', ln_trabbc
170         WRITE(numout,*) '      type of geothermal flux                      nn_geoflx     = ', nn_geoflx
171         WRITE(numout,*) '      Constant geothermal flux value               rn_geoflx_cst = ', rn_geoflx_cst
172         WRITE(numout,*)
173      ENDIF
174
175      IF( ln_trabbc ) THEN             !==  geothermal heating  ==!
176         !
177         ALLOCATE( qgh_trd0(jpi,jpj) )    ! allocation
178         !
179         SELECT CASE ( nn_geoflx )        ! geothermal heat flux / (rauO * Cp)
180         !
181         CASE ( 1 )                          !* constant flux
182            IF(lwp) WRITE(numout,*) '      *** constant heat flux  =   ', rn_geoflx_cst
183            qgh_trd0(:,:) = r1_rau0_rcp * rn_geoflx_cst
184            !
185         CASE ( 2 )                          !* variable geothermal heat flux : read the geothermal fluxes in mW/m2
186            IF(lwp) WRITE(numout,*) '      *** variable geothermal heat flux'
187            !
188            ALLOCATE( sf_qgh(1), STAT=ierror )
189            IF( ierror > 0 ) THEN
190               CALL ctl_stop( 'tra_bbc_init: unable to allocate sf_qgh structure' )   ;
191               RETURN
192            ENDIF
193            ALLOCATE( sf_qgh(1)%fnow(jpi,jpj,1)   )
194            IF( sn_qgh%ln_tint )ALLOCATE( sf_qgh(1)%fdta(jpi,jpj,1,2) )
195            ! fill sf_chl with sn_chl and control print
196            CALL fld_fill( sf_qgh, (/ sn_qgh /), cn_dir, 'tra_bbc_init',   &
197               &          'bottom temperature boundary condition', 'nambbc' )
198
199            CALL fld_read( nit000, 1, sf_qgh )                         ! Read qgh data
200            qgh_trd0(:,:) = r1_rau0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2
201            !
202         CASE DEFAULT
203            WRITE(ctmp1,*) '     bad flag value for nn_geoflx = ', nn_geoflx
204            CALL ctl_stop( ctmp1 )
205            !
206         END SELECT
207         !
208      ELSE
209         IF(lwp) WRITE(numout,*) '      *** no geothermal heat flux'
210      ENDIF
211      !
212   END SUBROUTINE tra_bbc_init
213
214   !!======================================================================
215END MODULE trabbc
Note: See TracBrowser for help on using the repository browser.