source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

File size: 5.4 KB
Line 
1MODULE insitu_tem
2
3  USE dom_oce  ! ocean space and time domain
4  USE oce,            ONLY: tsn
5  USE par_oce,        ONLY: jpi, jpj, jpk, jpkm1
6  USE lbclnk,         ONLY: lbc_lnk
7  USE lib_mpp         ! MPP library
8
9  USE yomhook, ONLY: lhook, dr_hook
10  USE parkind1, ONLY: jprb, jpim
11
12  IMPLICIT NONE
13  PRIVATE
14
15  PUBLIC theta2t
16
17  !! * Accessibility
18  PUBLIC insitu_tem_alloc          ! routines called by step.F90
19
20  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: insitu_t
21
22   !! * Substitutions
23#  include "domzgr_substitute.h90"
24
25CONTAINS
26
27   REAL FUNCTION insitu_tem_alloc()
28      !!----------------------------------------------------------------------
29      INTEGER, DIMENSION(2) :: ierr
30      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
31      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
32      REAL(KIND=jprb)               :: zhook_handle
33
34      CHARACTER(LEN=*), PARAMETER :: RoutineName='INSITU_TEM_ALLOC'
35
36      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
37
38      !!----------------------------------------------------------------------
39      !
40      ierr = 0
41      !
42      ALLOCATE( insitu_t(jpi,jpj,jpk), STAT=ierr(1) )
43         !
44      insitu_tem_alloc = MAXVAL(ierr)
45      IF( lk_mpp )   CALL mpp_sum( insitu_tem_alloc )
46      !
47      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
48  END FUNCTION insitu_tem_alloc
49
50!-----------------------------------------------------------------------------------
51!
52! Calculate the insitu temperature by integrating the adiabatic lapse rate from zero
53! to pressure at depth of tracer level. Based on UM subroutine POTTEM and function ATG
54!
55! Initial version: D. Acreman June 2006
56!
57!-----------------------------------------------------------------------------------
58
59  SUBROUTINE theta2t()
60
61    INTEGER, PARAMETER :: num_steps=10                ! number of steps in integration
62    INTEGER            :: step                        ! iteration counter
63    INTEGER            :: ji, jj, jk                  ! loop indices
64    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zP            ! pressure (decibars)
65    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zT            ! temperature at pressure p
66    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTB           ! temperature at p-dp
67    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTA           ! temperature at p+dp
68    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zDP           ! pressure step
69    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zSS           ! salinity(PSU) - 35.0
70    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zLAPSE        ! adiabatic lapse rate
71    INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
72    INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
73    REAL(KIND=jprb)               :: zhook_handle
74
75    CHARACTER(LEN=*), PARAMETER :: RoutineName='THETA2T'
76
77    IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
78
79
80!CDIR IEXPAND (ATG)
81
82! Set the pressure interval for the integration. The integration is carried out from
83! zero (pressure at the surface) to the pressure at the depth of the tracer point. The
84! pressure in decibars is represented by the depth in metres. Pressures are "Oceanographic"
85! pressures equal to absolute pressure minus one atmosphere
86     zDP(:,:,:) = 0.0
87     DO jk = 1, jpkm1
88        DO jj = 1, jpj
89           DO ji = 1, jpi
90              ! These loops expanded for case where fsdept may be 1D
91              zDP(ji,jj,jk) = fsdept(ji,jj,jk) / real(num_steps)
92           END DO
93        END DO
94     END DO
95
96! Salinity at each point
97     zSS(:,:,:)      = tsn(:,:,:,jp_sal) - 35.0 
98
99! Set initial values of temperature and pressure. zT is the temperature at pressure zP,
100! zTB is the temperature at pressure zP-zdP and zTA is the temperature at pressure zP+zdP
101     zT(:,:,:)  = tsn(:,:,:,jp_tem)
102     zP(:,:,:)  = 0.0         ! Pressure at surface
103     CALL ATG(zP, zT, zSS, zLAPSE)
104     zTB(:,:,:) = zT(:,:,:)  - zLAPSE(:,:,:) * zDP(:,:,:)
105
106     interation: DO step=1, num_steps
107        ! Calculate lapse rate (dT/dP) and hence TA
108        CALL ATG(zP, zT, zSS, zLAPSE)
109        zTA(:,:,:) = zTB(:,:,:) + 2.0 * zLAPSE(:,:,:) * zDP(:,:,:)
110        ! Have calculated TB, T and TA for this pressure, now advance solution by dP
111        zP(:,:,:)  = zP(:,:,:)  + zDP
112        zTB(:,:,:) = zT(:,:,:)
113        zT(:,:,:)  = zTA(:,:,:)
114
115     END DO interation
116
117     insitu_t(:,:,:) = zT(:,:,:) * tmask(:,:,:)
118     CALL lbc_lnk( insitu_t,  'T', 1.0)
119
120    IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
121   END SUBROUTINE theta2t
122
123   SUBROUTINE ATG(P,T,DS,zLAPSE)
124
125     REAL, INTENT(IN   ) :: P(jpi,jpj,jpk)        ! PRESSURE (DECIBARS)
126     REAL, INTENT(IN   ) :: T(jpi,jpj,jpk)        ! TEMPERATURE (DEG C)
127     REAL, INTENT(IN   ) :: DS(jpi,jpj,jpk)       ! SALINITY (PSU) -35.0
128     REAL, INTENT(  OUT) :: zLAPSE(jpi,jpj,jpk)   ! LAPSE RATE
129     INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
130     INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
131     REAL(KIND=jprb)               :: zhook_handle
132
133     CHARACTER(LEN=*), PARAMETER :: RoutineName='ATG'
134
135     IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
136
137
138       zLAPSE = ((( -2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P &
139           + ((2.7759E-12*T-1.1351E-10)*DS+((-5.4481E-14*T    &
140           + 8.733E-12)*T-6.7795E-10)*T+1.8741E-8))*P         &
141           + (-4.2393E-8*T+1.8932E-6)*DS                      &
142           + ((6.6228E-10*T-6.836E-8)*T+8.5258E-6)*T+3.5803E-5 
143
144     IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
145   END SUBROUTINE ATG
146
147END MODULE insitu_tem
Note: See TracBrowser for help on using the repository browser.