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.
paresp.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/paresp.F90 @ 3612

Last change on this file since 3612 was 3611, checked in by pabouttier, 11 years ago

Add TAM code and ORCA2_TAM configuration - see Ticket #1007
File size: 10.1 KB
Line 
1MODULE paresp
2   !!======================================================================
3   !!                       ***  MODULE paresp ***
4   !! NEMOVAR : Weights for an energy-type scalar product
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!  par_esp : Compute and store weights for an energy-type scalar product
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE par_kind
12   USE par_oce
13   USE eosbn2
14   USE phycst
15   USE zdf_oce
16   USE dom_oce
17   USE lib_mpp          ! MPP library
18   USE in_out_manager   ! I/O stuff
19   USE wrk_nemo         ! Memory Allocation
20
21   IMPLICIT NONE
22
23   !! * Routine accessibility
24   PRIVATE
25
26   PUBLIC &
27      & par_esp       !: Compute and store energy weights
28
29   !! * Share module variables
30   REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC ::   &
31      & wesp_t,   &   !: Normalized energy weights for temperature
32      & wesp_s        !: Normalized energy weights for salinity
33   REAL(wp), PUBLIC ::   &
34      & wesp_u,   &   !: Normalized energy weights for velocity
35      & wesp_ssh, &   !: Normalized energy weights for SSH
36      & wesp_tau, &   !: Normalized energy weights for windstress
37      & wesp_q,   &   !: Normalized energy weights for heat flux
38      & wesp_emp      !: Normalized energy weights for EmP
39
40CONTAINS
41   INTEGER FUNCTION par_esp_alloc()
42      !!----------------------------------------------------------------------
43      !!                ***  FUNCTION exa_mpl_alloc  ***
44      !!----------------------------------------------------------------------
45      ALLOCATE( wesp_t(jpk) , wesp_s(jpk) , STAT= par_esp_alloc )
46      !
47      IF( lk_mpp             )   CALL mpp_sum ( par_esp_alloc )
48      IF( par_esp_alloc /= 0 )   CALL ctl_warn('par_esp_alloc: failed to allocate arrays')
49      !
50   END FUNCTION par_esp_alloc
51
52   SUBROUTINE par_esp
53      !!-----------------------------------------------------------------------
54      !!
55      !!                  ***  ROUTINE par_esp  ***
56      !!
57      !! ** Purpose : Compute and store weights for an energy-type scalar
58      !!              product.
59      !!
60      !! ** Method  :
61      !!      The weighting coefficients are computed from an analytical
62      !!      linear density profile which is similar to the one used in
63      !!      default option in OPA direct.
64      !!
65      !!      The weighting coefficients for the forcing fields are estimated
66      !!      from the surface boundary conditions.
67      !!
68      !!      The weighting coefficients for the velocity variables are
69      !!      normalized to one.
70      !!
71      !!       wesp_t(k) = ( g / N(k) )^2 x alpha^2
72      !!       wesp_s(k) = ( g / N(k) )^2 x beta^2
73      !!       wesp_u    = 1
74      !!       wesp_ssh  = g
75      !!       wesp_tau  = ( ( e3w(1) / ( avu x rho0 ) )^2
76      !!       wesp_q    = ( ( g / N(1) x alpha
77      !!                  x e3w(1) / ( avT x rho0 x cp ) )^2
78      !!       wesp_emp  = ( ( g / N(1) x beta
79      !!                  x e3w(1) x S(1) / avT )^2
80      !!
81      !!      where
82      !!
83      !!        N(k)   is the Brunt-Vaisala frequency (depth-dependent)
84      !!        alpha  is the thermal expansion coefficient
85      !!        beta   is the salinity expansion coefficient
86      !!        g      is gravity
87      !!        avu    is vertical eddy viscosity coefficent for dynamics
88      !!        avT    is vertical eddy diffusivity coefficent for tracers
89      !!        e3w(1) is the vertical scale factor at top w-point
90      !!        rho0   is a reference density
91      !!        S(1)   is a reference salinity at the surface
92      !!
93      !!      The complete scalar product < ; > for the state vector
94      !!
95      !!       dx = ( du , dv , dT , dS , deta , dq , demp , dtaux , dtauy )
96      !!
97      !!      involves the volume/area elements:
98      !!
99      !!      < dx ; dx > = sum_i,j (  deta^2  * e1T * e2T           * wesp_ssh
100      !!                             + dq^2    * e1T * e2T * e3w(1)  * wesp_q
101      !!                             + demp^2  * e1T * e2T * e3w(1)  * wesp_emp
102      !!                             + dtaux^2 * e1u * e2u * e3uw(1) * wesp_tau
103      !!                             + dtauy^2 * e1v * e2v * e3vw(1) * wesp_tau
104      !!                    + sum_k (  du^2 * e1u * e2u * e3u * wesp_u
105      !!                             + dv^2 * e1v * e2v * e3v * wesp_u
106      !!                             + dT^2 * e1T * e2T * e3T * wesp_T
107      !!                             + dS^2 * e1T * e2T * e3T * wesp_S )    )
108      !!
109      !! ** Action  :
110      !!
111      !! History :
112      !!        ! 97-06 (A. Weaver, J. Vialard) OPAVAR version
113      !!        ! 07-11 (A. Weaver) NEMOVAR version based on OPAVAR paresp.F
114      !!        ! 2011-07 (J. While) Adapted to deal with 2D grids
115      !!-----------------------------------------------------------------------
116      !! * Modules used
117
118      !! * Arguments
119
120      !! * Local declarations
121      REAL(wp) :: zgrau0,  zk, zt0, zs0, zrau0
122      REAL(wp), POINTER, DIMENSION(:) :: ztan,  zsan
123      REAL(wp), POINTER, DIMENSION(:) :: zrhan, zbn2an
124      INTEGER :: jk, ierr
125      !!----------------------------------------------------------------------
126      !
127      ierr = par_esp_alloc()
128      CALL wrk_alloc( jpk, ztan, zsan, zrhan, zbn2an )
129      !
130      !--------------------------------------------------------------------
131      ! Local constant initialization
132      !--------------------------------------------------------------------
133
134      zt0   = 19.0_wp    ! Reference temperature
135      zs0   = 35.0_wp    ! Reference salinity
136      zrau0 = 1025.0_wp  ! Reference density
137
138      zgrau0 = -1.0_wp * grav / zrau0
139
140      !--------------------------------------------------------------------
141      ! Initialize rho with an analytical profile
142      !--------------------------------------------------------------------
143
144      !  Define analytical temperature profile
145
146      DO jk = 1, jpk
147         ztan(jk) =  7.5_wp * ( 1.0_wp - TANH( ( gdept_0(jk) - 80.0_wp ) &
148            &                 / 30.0_wp ) ) &
149            &      + 10.0_wp * ( 6000.0_wp - gdept_0(jk) ) / 6000.0_wp
150         ztan(jk) = ABS( ztan(jk) )
151      END DO
152#if defined key_pomme_r025
153      ztan(44) = 0.20
154      ztan(45) = 0.15
155      ztan(46) = 0.10
156#endif
157
158      !  Define analytical salinity profile
159
160      DO jk = 1, jpk
161         zsan(jk) = 35.50_wp
162      END DO
163
164      !  Define analytical density profile
165
166      DO jk = 1, jpk
167         zrhan(jk) =  zrau0 * ( 1.0_wp - rn_alpha * ( ztan(jk) - zt0 ) &
168            &                          + rn_beta  * ( zsan(jk) - zs0 ) )
169      END DO
170
171      !--------------------------------------------------------------------
172      ! Calculate N^2 at T and S points
173      !--------------------------------------------------------------------
174
175      IF ( jpk > 2 ) THEN
176         DO jk = 2, jpkm1
177            zbn2an(jk) = ( zgrau0 / 2.0_wp ) &
178               & * (   ( ( zrhan(jk-1) - zrhan(jk  ) ) / e3w_0(jk  ) ) &
179               &     + ( ( zrhan(jk  ) - zrhan(jk+1) ) / e3w_0(jk+1) ) )
180         END DO
181      ELSE
182         zbn2an(2) = 1._wp
183      ENDIF
184
185      zbn2an(1)   = zbn2an(2)
186      zbn2an(jpk) = zbn2an(jpkm1)
187
188      !--------------------------------------------------------------------
189      ! Calculate energy-type weights
190      !--------------------------------------------------------------------
191
192      DO jk = 1, jpk
193         zk = ( grav * grav ) / zbn2an(jk)
194         wesp_t(jk) = zk * rn_alpha * rn_alpha
195         wesp_s(jk) = zk * rn_beta  * rn_beta
196      END DO
197
198      wesp_u   = 1.0_wp
199      wesp_ssh = grav
200      wesp_tau = ( e3w_0(1) / ( rn_avm0 * zrau0 ) )**2
201      wesp_q   = wesp_t(1) * ( e3w_0(1) / ( rn_avt0 * zrau0 * rcp ) )**2
202      wesp_emp = wesp_s(1) * ( e3w_0(1) * zsan(1) / rn_avt0  )**2
203
204      !--------------------------------------------------------------------
205      ! Print
206      !--------------------------------------------------------------------
207
208      IF (lwp) THEN
209         WRITE(numout,*)
210         WRITE(numout,*) ' par_esp: Analytical T, S, rho, N^2 profiles'
211         WRITE(numout,*) ' -------'
212         WRITE(numout,*)
213         WRITE(numout,995)
214         WRITE(numout,996)
215         DO jk = 1, jpk
216            IF(lwp)WRITE(numout,1007) jk, ztan(jk), zsan(jk), &
217               &                      zrhan(jk), zbn2an(jk)
218         END DO
219         WRITE(numout,*)
220         WRITE(numout,*) ' par_esp: Energy scalar product weights', &
221            &            ' for T and S'
222         WRITE(numout,*) ' -------'
223         WRITE(numout,*)
224         WRITE(numout,998)
225         WRITE(numout,999)
226         DO jk = 1, jpk
227            WRITE(numout,1006) jk, wesp_t(jk), wesp_s(jk)
228         END DO
229         WRITE(numout,*)
230         WRITE(numout,*) ' par_esp: Energy scalar product weights', &
231            &            ' for u, SSH, tau, Q and EmP'
232         WRITE(numout,*) ' -------'
233         WRITE(numout,*)
234         WRITE(numout,*) '          wesp_u   = ', wesp_u
235         WRITE(numout,*) '          wesp_ssh = ', wesp_ssh
236         WRITE(numout,*) '          wesp_tau = ', wesp_tau
237         WRITE(numout,*) '          wesp_q   = ', wesp_q
238         WRITE(numout,*) '          wesp_emp = ', wesp_emp
239
240 995     FORMAT('    level        T              S    ', &
241            &   '        rho             N^2 ')
242 996     FORMAT('    -----     --------       --------', &
243            &   '     ----------     ----------')
244 998     FORMAT('    level   wesp_T        wesp_S   ')
245 999     FORMAT('    ----- ----------   ------------')
246 1005    FORMAT(4X,I2,3X,E12.5,1X,E12.5)
247 1006    FORMAT(4X,I2,3X,E12.5,1X,E12.5)
248 1007    FORMAT(4X,I2,6X,F10.5,5X,F10.5,5X,F10.5,5X,E12.5)
249         CALL flush( numout )
250
251      ENDIF
252
253      DO jk = 1, jpk
254         IF ( zbn2an(jk) < 0.0_wp ) THEN
255            IF(lwp)WRITE(numout,990) zbn2an(jk), jk
256990         FORMAT(' paresp : Error unstable density profile;', &
257               &   ' zbn2an = ',E12.5,' at level ',I3)
258!            CALL ctl_stop( 'paresp; unstable density profile' )
259             IF(lwp)WRITE(numout,*)'WARNING'
260        IF(lwp)WRITE(numout,*)'paresp; unstable density profile'
261         ENDIF
262      ENDDO
263      CALL wrk_dealloc( jpk, ztan, zsan, zrhan, zbn2an )
264
265   END SUBROUTINE par_esp
266
267END MODULE paresp
Note: See TracBrowser for help on using the repository browser.