source: NEMO/trunk/src/OCE/DIU/diurnal_bulk.F90 @ 10069

Last change on this file since 10069 was 10069, checked in by nicolasmartin, 2 years ago

Fix mistakes of previous commit on SVN keywords property

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