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.
diurnal_bulk.F90 in NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/DIU – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/DIU/diurnal_bulk.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 11.1 KB
Line 
1MODULE diurnal_bulk
2   !!======================================================================
3   !!                    ***  MODULE  diurnal_bulk  ***
4   !!     Takaya model of diurnal warming (Takaya, 2010)
5   !!=====================================================================
6   !! History :        !  11-10  (J. While)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   diurnal_sst_bulk_init  : initialise diurnal model
11   !!   diurnal_sst_bulk_step  : time-step the diurnal model
12   !!----------------------------------------------------------------------
13   USE par_kind
14   USE phycst
15   USE dom_oce
16   USE lib_mpp
17   USE solfrac_mod
18   USE in_out_manager
19   
20   IMPLICIT NONE
21   PRIVATE
22   
23   ! Namelist parameters
24   LOGICAL, PUBLIC :: ln_diurnal
25   LOGICAL, PUBLIC :: ln_diurnal_only
26
27   ! Parameters
28   REAL(wp), PRIVATE, PARAMETER :: pp_alpha = 2.0e-4_wp
29   REAL(wp), PRIVATE, PARAMETER :: pp_veltol = 0._wp
30   REAL(wp), PRIVATE, PARAMETER :: pp_min_fvel = 1.e-10_wp 
31   
32   ! Key variables
33   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_dsst     ! Delta SST
34   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_solfrac  ! Fraction of
35   !                                                           ! absorbed radiation
36
37   PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step
38   
39   !!----------------------------------------------------------------------
40CONTAINS
41   
42   SUBROUTINE diurnal_sst_bulk_init
43      !!----------------------------------------------------------------------
44      !! *** ROUTINE diurnal_sst_init ***
45      !!
46      !! ** Purpose : Initialise the Takaya diurnal model   
47      !!----------------------------------------------------------------------     
48      INTEGER ::   ios   ! local integer
49      !!
50      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
51      !!----------------------------------------------------------------------     
52
53      ! Read the namelist
54      READ  ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 )
55901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdiu in reference namelist' )
56      READ  ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
57902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdiu in configuration namelist' )     
58      !
59      IF( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN
60         CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" )
61      ENDIF
62     
63      IF( ln_diurnal ) THEN     
64         !         
65         ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
66         !
67         x_solfrac = 0._wp         ! Initialise the solar fraction
68         x_dsst    = 0._wp
69         !
70         IF( ln_diurnal_only ) THEN
71            CALL ctl_warn( "ln_diurnal_only set; only the diurnal component of SST will be calculated" )
72         ENDIF
73      ENDIF
74     
75   END SUBROUTINE diurnal_sst_bulk_init
76
77
78   SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, p_rdt,   &
79            &                  pla, pthick, pcoolthick, pmu, &
80            &                  p_fvel_bkginc, p_hflux_bkginc)
81      !!----------------------------------------------------------------------
82      !! *** ROUTINE diurnal_sst_takaya_step ***
83      !!
84      !! ** Purpose :   Time-step the Takaya diurnal model
85      !!
86      !! ** Method :    1) Calculate the Obukhov length
87      !!                2) Calculate the Similarity function
88      !!                2) Calculate the increment to dsst
89      !!                3) Apply the increment
90      !! ** Reference : Refinements to a prognostic scheme of skin sea surface
91      !!                temperature, Takaya et al, JGR, 2010
92      !!----------------------------------------------------------------------
93      INTEGER                               , INTENT(in) ::   kt             ! time step
94      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   psolflux       ! solar flux (Watts)
95      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   pqflux         ! heat (non-solar) flux (Watts)
96      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   ptauflux       ! wind stress  (kg/ m s^2)
97      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   prho           ! water density  (kg/m^3)
98      REAL(wp)                              , INTENT(in) ::   p_rdt          ! time-step
99      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pLa            ! Langmuir number
100      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pthick         ! warm layer thickness (m)
101      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pcoolthick     ! cool skin thickness (m)
102      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pmu            ! mu parameter
103      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   p_hflux_bkginc ! increment to the heat flux
104      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   p_fvel_bkginc  ! increment to the friction velocity
105      !
106      INTEGER :: ji,jj
107      LOGICAL  :: ll_calcfrac
108      REAL(wp), DIMENSION(jpi,jpj) :: z_fvel              ! friction velocity     
109      REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
110      REAL(wp), DIMENSION(jpi,jpj) :: z_abflux            ! absorbed flux           
111      REAL(wp), DIMENSION(jpi,jpj) :: z_fla               ! Langmuir function value
112      !!----------------------------------------------------------------------
113
114      ! Set optional arguments to their defaults
115      IF( .NOT. PRESENT( pthick )   ) THEN   ;   zthick(:,:) = 3._wp
116      ELSE                                   ;   zthick(:,:) = pthick(:,:)
117      ENDIF
118      IF( .NOT. PRESENT(pcoolthick) ) THEN   ;   zcoolthick(:,:) = 0._wp
119      ELSE                                   ;   zcoolthick(:,:) = pcoolthick(:,:)
120      ENDIF
121      IF( .NOT. PRESENT( pmu )      ) THEN   ;   zmu(:,:) = 0.3_wp
122      ELSE                                   ;   zmu(:,:) = pmu(:,:)
123      ENDIF
124      IF( .NOT. PRESENT(pla) ) THEN          ;   zla(:,:) = 0.3_wp
125      ELSE                                   ;   zla(:,:) = pla(:,:)
126      ENDIF
127     
128      ! If not done already, calculate the solar fraction
129      IF ( kt==nit000 ) THEN
130         DO jj = 1,jpj
131            DO ji = 1, jpi
132               IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
133                  &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
134            END DO
135         END DO   
136      ENDIF   
137
138      ! convert solar flux and heat flux to absorbed flux   
139      WHERE ( tmask(:,:,1) == 1._wp) 
140         z_abflux(:,:) = (  x_solfrac(:,:) *  psolflux (:,:)) + pqflux(:,:)     
141      ELSEWHERE
142         z_abflux(:,:) = 0._wp
143      ENDWHERE
144      IF( PRESENT(p_hflux_bkginc) ) z_abflux(:,:) = z_abflux(:,:) + p_hflux_bkginc   ! Optional increment
145      WHERE ( ABS( z_abflux(:,:) ) < rsmall )
146         z_abflux(:,:) = rsmall
147      ENDWHERE 
148     
149      ! Calculate the friction velocity
150      WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )   
151         z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
152      ELSEWHERE
153         z_fvel(:,:) = 0._wp
154      ENDWHERE
155      IF( PRESENT(p_fvel_bkginc) ) z_fvel(:,:) = z_fvel(:,:) + p_fvel_bkginc   ! Optional increment
156     
157       
158       
159      ! Calculate the Langmuir function value
160      WHERE ( tmask(:,:,1) == 1.)
161         z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) ) 
162      ELSEWHERE
163         z_fla(:,:) = 0._wp
164      ENDWHERE     
165     
166      ! Increment the temperature using the implicit solution
167      x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:),   &
168         &                       z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:) )
169      !
170   END SUBROUTINE diurnal_sst_takaya_step
171
172   
173   FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, &
174                          p_fla, pmu, pthick, prho )
175                         
176      IMPLICIT NONE
177     
178      ! Function definition
179      REAL(wp), DIMENSION(jpi,jpj) :: t_imp
180      ! Dummy variables
181      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst     ! Delta SST
182      REAL(wp), INTENT(IN)                     :: p_rdt      ! Time-step
183      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux   ! Heat forcing
184      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel     ! Friction velocity
185      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla      ! Langmuir number
186      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu        ! Structure parameter
187      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick     ! Layer thickness
188      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho       ! Water density
189   
190      ! Local variables
191      REAL(wp) :: z_olength          ! Obukhov length
192      REAL(wp) :: z_sigma, z_sigma2
193      REAL(wp) :: z_term1, z_term2     
194      REAL(wp) :: z_stabfunc         ! stability function value
195      REAL(wp) :: z_fvel     
196     
197      CHARACTER(200) :: warn_string
198       
199      INTEGER :: ji,jj
200                     
201      DO jj = 1, jpj
202         DO ji = 1, jpi   
203           
204            ! Only calculate outside tmask
205            IF ( tmask(ji,jj,1) /= 1._wp ) THEN
206               t_imp(ji,jj) = 0._wp
207               CYCLE   
208            END IF
209           
210            IF (p_fvel(ji,jj) < pp_min_fvel) THEN
211               z_fvel = pp_min_fvel
212               WRITE(warn_string,*) "diurnal_sst_takaya step: "&
213               &//"friction velocity < minimum\n" &
214               &//"Setting friction velocity =",pp_min_fvel
215               CALL ctl_warn(warn_string)
216               
217            ELSE
218               z_fvel = p_fvel(ji,jj)
219            ENDIF
220                 
221            ! Calculate the Obukhov length
222            IF ( (z_fvel < pp_veltol ) .AND. &
223            &    (p_dsst(ji,jj) > 0._wp) ) THEN
224               z_olength =  z_fvel  / &
225               &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
226               &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
227            ELSE
228               z_olength = &
229               &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
230               &   ( vkarmn * grav * pp_alpha *&
231               &     p_abflux(ji,jj) )
232            ENDIF
233             
234            ! Calculate the stability function         
235            z_sigma = pthick(ji,jj) / z_olength
236            z_sigma2 = z_sigma * z_sigma
237 
238            IF ( z_sigma >= 0. ) THEN
239               z_stabfunc = 1._wp + &
240               &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
241               &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
242               &    z_sigma2 ) )
243            ELSE
244               z_stabfunc = 1._wp / &
245               &                   SQRT( 1._wp - 16._wp * z_sigma )
246            ENDIF
247
248            ! Calculate the T increment
249            z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
250            & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
251           
252             
253            z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
254            &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
255            &      ( pthick(ji,jj) * z_stabfunc ) )     
256         
257            t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
258                           ( 1._wp - p_rdt * z_term2 )
259
260         END DO
261      END DO
262     
263      END FUNCTION t_imp
264
265END MODULE diurnal_bulk
Note: See TracBrowser for help on using the repository browser.