source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90 @ 6010

Last change on this file since 6010 was 6010, checked in by timgraham, 5 years ago

Tidying of DIU code

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