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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DIU – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DIU/diu_bulk.F90 @ 11463

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

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert NST routines in preparation for getting AGRIF back up and running. AGRIF conv stage now works but requires some renaming of recently changes DIU modules (included in this commit). AGRIF compile and link stage not yet working (agrif routines need to be passed the time-level indices) but non-AGRIF SETTE tests are all OK
  • Property svn:keywords set to Id
File size: 11.1 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
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      REWIND( numnam_ref )
55      READ  ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 )
56901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdiu in reference namelist', lwp )
57      REWIND( numnam_cfg )
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
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 jj = 1,jpj
133            DO ji = 1, jpi
134               IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
135                  &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
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
169      x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:),   &
170         &                       z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:) )
171      !
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 diu_bulk
Note: See TracBrowser for help on using the repository browser.