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 branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU – NEMO

source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90 @ 11137

Last change on this file since 11137 was 11137, checked in by jcastill, 5 years ago

Add missing files
File size: 12.0 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  : initalise diurnal model
11   !!   diurnal_sst_bulk_step  : timestep 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   
22   !Namelist parameters
23   LOGICAL, PUBLIC :: ln_diurnal = .FALSE.
24   LOGICAL, PUBLIC :: ln_diurnal_only = .FALSE.
25
26   !Parameters
27   REAL(wp), PRIVATE, PARAMETER :: pp_alpha = 2.0e-4_wp
28   REAL(wp), PRIVATE, PARAMETER :: pp_veltol = 0._wp
29   REAL(wp), PRIVATE, PARAMETER :: pp_min_fvel = 1.e-10_wp 
30   
31   !Key variables
32   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_dsst     !Delta SST
33   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_solfrac  !Fraction of
34                                                      !absorbed
35                                                      !radiation
36
37   PRIVATE
38   PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step
39   
40   CONTAINS
41   
42   SUBROUTINE diurnal_sst_bulk_init
43      !!----------------------------------------------------------------------
44      !! *** ROUTINE diurnal_sst_init ***
45      !!
46      !! ** Purpose : Initalise the Takaya diurnal model
47     
48      !!----------------------------------------------------------------------
49     
50      IMPLICIT NONE
51
52      INTEGER ::   ios                 ! Local integer output status for namelist read
53     
54      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
55     
56      !Read the namelist
57      REWIND ( numnam_ref )              ! Read Namelist namdiu in reference namelist
58      READ   ( numnam_ref, namdiu, IOSTAT=ios, ERR= 901 )
59901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in reference namelist', lwp )
60
61      REWIND( numnam_cfg )              ! Namelist namdiu in configuration namelist
62      READ  ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
63902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in configuration namelist', lwp )
64      IF(lwm) WRITE ( numond, namdiu )
65     
66      IF(lwp) THEN                   ! Control print
67         WRITE(numout,*)
68         WRITE(numout,*) 'namdiu_init : Diurnal Model settings are:'
69         WRITE(numout,*)    'ln_diurnal = ', ln_diurnal
70         WRITE(numout,*)    'ln_diurnal_only = ', ln_diurnal_only
71      ENDIF
72
73      IF ( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN
74         CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" )
75      ENDIF
76     
77      IF ( ln_diurnal ) THEN     
78         
79         !Allocate arrays
80         ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
81         
82         !Initalise the solar fraction
83         x_solfrac=0._wp
84
85         IF ( ln_diurnal_only ) THEN
86            CALL ctl_warn( "ln_diurnal_only set; only the diurnal component of SST will be calculated" )
87         ENDIF
88      ENDIF
89     
90   END SUBROUTINE diurnal_sst_bulk_init
91   
92   SUBROUTINE diurnal_sst_takaya_step(psolflux, pqflux, ptauflux, prho, p_rdt,&
93            &                  pla, pthick, pcoolthick, pmu, ld_calcfrac, &
94            &                  p_fvel_bkginc, p_hflux_bkginc)
95      !!----------------------------------------------------------------------
96      !! *** ROUTINE diurnal_sst_takaya_step ***
97      !!
98      !! ** Purpose :   Timestep the Takaya diurnal model
99      !!
100      !! ** Method :    1) Calculate the Obukhov length
101      !!                2) Calculate the Simularity function
102      !!                2) Calculate the increment to dsst
103      !!                3) Apply the increment
104      !! ** Reference : Refinements to a prognostic scheme of skin sea surface
105      !!                temperature, Takaya et al, JGR, 2010
106      !!----------------------------------------------------------------------
107     
108      IMPLICIT NONE
109     
110      !Dummy variables
111      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psolflux  !solar flux (Watts)
112      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pqflux    !heat (non-solar)
113                                                            !flux (Watts)
114      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: ptauflux  !wind stress
115                                                            !(kg/ m s^2)
116      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: prho      !water density
117                                                            !(kg/m^3)
118      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pLa   
119                                                     !Langmuir number
120      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pthick 
121                                                     !warm layer thickness (m)
122      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pcoolthick 
123                                                     !cool skin thickness (m)
124      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pmu   
125                                                     !mu parameter
126      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: p_hflux_bkginc 
127                                                     !optional increment to the
128                                                     !heat flux
129      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: p_fvel_bkginc 
130                                                     !optional increment to the
131                                                     !friction velocity
132      REAL(wp), INTENT(IN) :: p_rdt              !timestep
133     
134      LOGICAL, OPTIONAL, &
135      &     INTENT(IN)  :: ld_calcfrac !Set TRUE to recalculate the
136                                               !solar fraction
137     
138      !Local variables
139      REAL(wp), DIMENSION(jpi,jpj) :: z_fvel              ! friction velocity     
140      REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
141      REAL(wp), DIMENSION(jpi,jpj) :: z_abflux            ! absorbed flux           
142      REAL(wp), DIMENSION(jpi,jpj) :: z_fla               !Langmuir fuction value
143     
144      LOGICAL  :: ll_calcfrac
145     
146      INTEGER :: ji,jj
147
148      !Set optional arguments to their defualts
149      IF ( .NOT. PRESENT(pthick) ) THEN
150         zthick(:,:) = 3._wp
151      ELSE
152         zthick(:,:) = pthick(:,:)
153      ENDIF
154      IF ( .NOT. PRESENT(pcoolthick) ) THEN
155         zcoolthick(:,:) = 0._wp
156      ELSE
157         zcoolthick(:,:) = pcoolthick(:,:)
158      ENDIF
159      IF ( .NOT. PRESENT(pmu) ) THEN
160         zmu(:,:) = 0.3_wp
161      ELSE
162         zmu(:,:) = pmu(:,:)
163      ENDIF
164      IF ( .NOT. PRESENT(pla) ) THEN
165         zla(:,:) = 0.3_wp
166      ELSE
167         zla(:,:) = pla(:,:)
168      ENDIF
169      IF ( .NOT. PRESENT(ld_calcfrac) ) THEN
170         ll_calcfrac = .FALSE.
171      ELSE
172         ll_calcfrac = ld_calcfrac
173      ENDIF     
174     
175      !If not done already, calculate the solar fraction
176      IF (ll_calcfrac ) THEN
177         DO jj = 1,jpj
178            DO ji = 1, jpi
179               IF(  ( x_solfrac(ji,jj) == 0. ) .AND. ( tmask(ji,jj,1) == 1. ) ) &
180               &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
181            END DO
182         END DO   
183      ENDIF   
184
185      !convert solar flux and heat flux to absorbed flux   
186      WHERE ( tmask(:,:,1) == 1.) 
187         z_abflux(:,:) = (  x_solfrac(:,:) *  psolflux (:,:)) + pqflux(:,:)     
188      ELSEWHERE
189         z_abflux(:,:) = 0._wp
190      ENDWHERE
191      IF( PRESENT(p_hflux_bkginc) ) z_abflux(:,:) = z_abflux(:,:) + p_hflux_bkginc !Optional increment
192      WHERE ( ABS( z_abflux(:,:) ) < rsmall )
193         z_abflux(:,:) = rsmall
194      ENDWHERE 
195     
196      !Calculate the friction velocity
197      WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )   
198         z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
199      ELSEWHERE
200         z_fvel(:,:) = 0.
201      ENDWHERE
202      IF( PRESENT(p_hflux_bkginc) ) z_fvel(:,:) = z_fvel(:,:) + p_fvel_bkginc !Optional increment
203     
204       
205       
206      !Calculate the Langmuir function value
207      WHERE ( tmask(:,:,1) == 1.)
208         z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) ) 
209      ENDWHERE     
210     
211      !!Increment the temperature using the implicit solution
212      x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:), &
213                                 z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:))
214     
215                               
216     
217   END SUBROUTINE diurnal_sst_takaya_step
218
219   
220   FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, &
221                          p_fla, pmu, pthick, prho )
222                         
223      IMPLICIT NONE
224     
225      !Function definition
226      REAL(wp), DIMENSION(jpi,jpj) :: t_imp
227      !Dummy variables
228      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst    !Delta SST
229      REAL(wp), INTENT(IN) ::                 p_rdt           !Timestep
230      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux  !Heat forcing
231      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel    !Friction velocity
232      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla     !Langmuir number
233      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu       !Structure parameter
234      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick    !Layer thickness
235      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho      !Water density
236   
237      !Local variables
238      REAL(wp) :: z_olength          ! Obukhov length
239      REAL(wp) :: z_sigma, z_sigma2
240      REAL(wp) :: z_term1, z_term2     
241      REAL(wp) :: z_stabfunc         ! stability function value
242      REAL(wp) :: z_fvel     
243     
244      CHARACTER(200) :: warn_string
245       
246      INTEGER :: ji,jj
247                     
248      DO jj = 1, jpj
249         DO ji = 1, jpi   
250           
251            !Only calculate outside tmask
252            IF ( tmask(ji,jj,1) /= 1. ) THEN
253               t_imp(ji,jj) = 0._wp
254               CYCLE   
255            END IF
256           
257            IF (p_fvel(ji,jj) < pp_min_fvel) THEN
258               z_fvel = pp_min_fvel
259               WRITE(warn_string,*) "diurnal_sst_takaya step: "&
260               &//"friction velocity < minimum\n" &
261               &//"Setting friction velocity =",pp_min_fvel
262               CALL ctl_warn(warn_string)
263               
264            ELSE
265               z_fvel = p_fvel(ji,jj)
266            ENDIF
267                 
268            !Calculate the Obukhov length
269            IF ( (z_fvel < pp_veltol ) .AND. &
270            &    (p_dsst(ji,jj) > 0.) ) THEN
271               z_olength =  z_fvel  / &
272               &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
273               &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
274            ELSE
275               z_olength = &
276               &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
277               &   ( vkarmn * grav * pp_alpha *&
278               &     p_abflux(ji,jj) )
279            ENDIF
280             
281            !Calculate the stability function         
282            z_sigma = pthick(ji,jj) / z_olength
283            z_sigma2 = z_sigma * z_sigma
284 
285            IF ( z_sigma >= 0 ) THEN
286               z_stabfunc = 1._wp + &
287               &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
288               &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
289               &    z_sigma2 ) )
290            ELSE
291               z_stabfunc = 1._wp / &
292               &                   SQRT( 1._wp - 16._wp * z_sigma )
293            ENDIF
294
295            !Calculate the T increment
296            z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
297            & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
298           
299             
300            z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
301            &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
302            &      ( pthick(ji,jj) * z_stabfunc ) )     
303         
304            t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
305                           ( 1._wp - p_rdt * z_term2 )
306
307         END DO
308      END DO
309     
310      END FUNCTION t_imp
311
312END MODULE diurnal_bulk
Note: See TracBrowser for help on using the repository browser.