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.
diu_bulk.F90 in NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DIU – NEMO

source: NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DIU/diu_bulk.F90 @ 13886

Last change on this file since 13886 was 13886, checked in by rlod, 3 years ago

phasing with trunk at revision r13787
  • Property svn:keywords set to Id
File size: 11.0 KB
Line 
1MODULE diu_bulk
2   !!======================================================================
3   !!                    ***  MODULE  diu_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      = .false.   ! force definition if diurnal_sst_bulk_init is not called
25   LOGICAL, PUBLIC :: ln_diurnal_only = .false.   ! force definition if diurnal_sst_bulk_init is not called
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   !! * Substitutions
39#  include "do_loop_substitute.h90"
40   
41   !!----------------------------------------------------------------------
42CONTAINS
43   
44   SUBROUTINE diurnal_sst_bulk_init
45      !!----------------------------------------------------------------------
46      !! *** ROUTINE diurnal_sst_init ***
47      !!
48      !! ** Purpose : Initialise the Takaya diurnal model   
49      !!----------------------------------------------------------------------     
50      INTEGER ::   ios   ! local integer
51      !!
52      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
53      !!----------------------------------------------------------------------     
54
55      ! Read the namelist
56      READ  ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 )
57901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdiu in reference namelist' )
58      READ  ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
59902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdiu in configuration namelist' )     
60      !
61      IF( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN
62         CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" )
63      ENDIF
64     
65      IF( ln_diurnal ) THEN     
66         !         
67         ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
68         !
69         x_solfrac = 0._wp         ! Initialise the solar fraction
70         x_dsst    = 0._wp
71         !
72         IF( ln_diurnal_only ) THEN
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
78
79
80   SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, p_rdt,   &
81            &                  pla, pthick, pcoolthick, pmu, &
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      !!----------------------------------------------------------------------
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
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
114      !!----------------------------------------------------------------------
115
116      ! Set optional arguments to their defaults
117      IF( .NOT. PRESENT( pthick )   ) THEN   ;   zthick(:,:) = 3._wp
118      ELSE                                   ;   zthick(:,:) = pthick(:,:)
119      ENDIF
120      IF( .NOT. PRESENT(pcoolthick) ) THEN   ;   zcoolthick(:,:) = 0._wp
121      ELSE                                   ;   zcoolthick(:,:) = pcoolthick(:,:)
122      ENDIF
123      IF( .NOT. PRESENT( pmu )      ) THEN   ;   zmu(:,:) = 0.3_wp
124      ELSE                                   ;   zmu(:,:) = pmu(:,:)
125      ENDIF
126      IF( .NOT. PRESENT(pla) ) THEN          ;   zla(:,:) = 0.3_wp
127      ELSE                                   ;   zla(:,:) = pla(:,:)
128      ENDIF
129     
130      ! If not done already, calculate the solar fraction
131      IF ( kt==nit000 ) THEN
132         DO_2D( 1, 1, 1, 1 )
133            IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
134               &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
135         END_2D
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_2D( 1, 1, 1, 1 )
202         
203         ! Only calculate outside tmask
204         IF ( tmask(ji,jj,1) /= 1._wp ) THEN
205            t_imp(ji,jj) = 0._wp
206            CYCLE   
207         END IF
208         
209         IF (p_fvel(ji,jj) < pp_min_fvel) THEN
210            z_fvel = pp_min_fvel
211            WRITE(warn_string,*) "diurnal_sst_takaya step: "&
212            &//"friction velocity < minimum\n" &
213            &//"Setting friction velocity =",pp_min_fvel
214            CALL ctl_warn(warn_string)
215           
216         ELSE
217            z_fvel = p_fvel(ji,jj)
218         ENDIF
219             
220         ! Calculate the Obukhov length
221         IF ( (z_fvel < pp_veltol ) .AND. &
222         &    (p_dsst(ji,jj) > 0._wp) ) THEN
223            z_olength =  z_fvel  / &
224            &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
225            &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
226         ELSE
227            z_olength = &
228            &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
229            &   ( vkarmn * grav * pp_alpha *&
230            &     p_abflux(ji,jj) )
231         ENDIF
232         
233         ! Calculate the stability function         
234         z_sigma = pthick(ji,jj) / z_olength
235         z_sigma2 = z_sigma * z_sigma
236         IF ( z_sigma >= 0. ) THEN
237            z_stabfunc = 1._wp + &
238            &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
239            &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
240            &    z_sigma2 ) )
241         ELSE
242            z_stabfunc = 1._wp / &
243            &                   SQRT( 1._wp - 16._wp * z_sigma )
244         ENDIF
245
246         ! Calculate the T increment
247         z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
248         & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
249         
250         
251         z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
252         &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
253         &      ( pthick(ji,jj) * z_stabfunc ) )     
254       
255         t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
256                        ( 1._wp - p_rdt * z_term2 )
257
258      END_2D
259     
260      END FUNCTION t_imp
261
262END MODULE diu_bulk
Note: See TracBrowser for help on using the repository browser.