Changeset 2188 for branches/dev_r2174_DCY/NEMO/OPA_SRC/SBC/sbcdcy.F90
- Timestamp:
- 2010-10-08T10:32:36+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_r2174_DCY/NEMO/OPA_SRC/SBC/sbcdcy.F90
r2187 r2188 4 4 !! Ocean forcing: compute the diurnal cycle 5 5 !!====================================================================== 6 !! History : 8.2! 2005-02 (D. Bernie) Original code7 !! 9.0 ! 2006-02 (S. Masson, G. Madec) adaptation to OPA98 !! 3.1 ! 2009-07 (J.M. Molines) adaptation to nemo3.16 !! History : OPA ! 2005-02 (D. Bernie) Original code 7 !! NEMO 2.0 ! 2006-02 (S. Masson, G. Madec) adaptation to NEMO 8 !! 3.1 ! 2009-07 (J.M. Molines) adaptation to v3.1 9 9 !!---------------------------------------------------------------------- 10 10 … … 20 20 IMPLICIT NONE 21 21 PRIVATE 22 INTEGER :: idayqsr ! day when parameters were computed 23 REAL(wp), DIMENSION(jpi,jpj) :: zaaa, zbbb, zccc, zab, ztmd, zdawn, zdusk, zscal ! parameters to compute the diurnal cycle 24 REAL(wp), DIMENSION(jpi,jpj) :: qsr_daily ! to hold daily mean QSR 22 INTEGER :: nday_qsr ! day when parameters were computed 23 REAL(wp), DIMENSION(jpi,jpj) :: raa , rbb , rcc , rab ! parameters used to compute the diurnal cycle 24 REAL(wp), DIMENSION(jpi,jpj) :: rtmd, rdawn, rdusk, rscal ! - - - - - 25 REAL(wp), DIMENSION(jpi,jpj) :: qsr_daily ! to hold daily mean QSR 25 26 26 PUBLIC sbc_dcy! routine called by sbc27 28 !!---------------------------------------------------------------------- 29 !! NEMO/OPA 3.3 , LOCEAN-IPSL(2010)27 PUBLIC sbc_dcy ! routine called by sbc 28 29 !!---------------------------------------------------------------------- 30 !! NEMO/OPA 3.3 , NEMO-consortium (2010) 30 31 !! $Id$ 31 32 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 32 33 !!---------------------------------------------------------------------- 33 34 34 CONTAINS 35 35 … … 40 40 !! ** Purpose : introduce a diurnal cycle of qsr from daily values 41 41 !! 42 !! ** Method : see Appendix A of 43 !! Bernie, DJ, Guilyardi, E, Madec, G, Slingo, JM and Woolnough, SJ 44 !! Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. Part 1: a diurnally forced OGCM 45 !! Climate Dynamics 29:6, 575-590 (2007) 42 !! ** Method : see Appendix A of Bernie et al. 2007. 46 43 !! 47 44 !! ** Action : redistribute daily QSR on each time step following the diurnal cycle 45 !! 46 !! reference : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007 47 !! Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. 48 !! Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590. 48 49 !!---------------------------------------------------------------------- 49 INTEGER, INTENT( in ) :: kt ! ocean time-step index 50 REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: pqsr ! QSR flux with diurnal cycle 51 !! 52 INTEGER :: ji, jj ! dummy loop indices 53 REAL(wp) :: fintegral, pt1, pt2, paaa, pbbb, pccc ! 50 INTEGER, INTENT(in ) :: kt ! ocean time-step index 51 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqsr ! QSR flux with diurnal cycle 52 !! 53 INTEGER :: ji, jj ! dummy loop indices 54 54 REAL(wp) :: ztwopi, zinvtwopi, zconvrad 55 55 REAL(wp) :: zlo, zup, zlousd, zupusd 56 56 REAL(wp) :: zdsws, zdecrad, ztx 57 57 REAL(wp) :: ztmp, ztmp1, ztmp2, ztest 58 !!--------------------------------------------------------------------- 59 60 !---------------------------------------------------------------------- 61 ! statement functions 62 63 fintegral(pt1, pt2, paaa, pbbb, pccc) = & 58 !---------------------------statement functions------------------------ 59 REAL(wp) :: fintegral, pt1, pt2, paaa, pbbb, pccc ! dummy statement function arguments 60 fintegral( pt1, pt2, paaa, pbbb, pccc ) = & 64 61 & paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2) & 65 62 & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1) 66 ! ----------------------------------------------------------------------63 !!--------------------------------------------------------------------- 67 64 68 65 ! Initialization … … 77 74 78 75 ! 79 IF (kt == nit000) THEN 80 ! 76 IF( kt == nit000 ) THEN ! first time step only 81 77 IF(lwp) THEN 82 78 WRITE(numout,*) … … 85 81 WRITE(numout,*) 86 82 ENDIF 87 idayqsr = 088 ! Compute Cneeded to compute the time integral of the diurnal cycle89 zccc(:,:) = zconvrad * glamt(:,:) - rpi83 nday_qsr = 0 84 ! Compute rcc needed to compute the time integral of the diurnal cycle 85 rcc(:,:) = zconvrad * glamt(:,:) - rpi 90 86 ! time of midday 91 ztmd(:,:) = 0.5 - glamt(:,:) / 360.92 ztmd(:,:) = MOD((ztmd(:,:) + 1.), 1.)87 rtmd(:,:) = 0.5 - glamt(:,:) / 360. 88 rtmd(:,:) = MOD( (rtmd(:,:) + 1.), 1. ) 93 89 ENDIF 94 90 … … 99 95 100 96 ! nday is the number of days since the beginning of the current month 101 IF( idayqsr /= nday ) THEN97 IF( nday_qsr /= nday ) THEN 102 98 ! save the day of the year and the daily mean of qsr 103 idayqsr = nday99 nday_qsr = nday 104 100 ! number of days since the previous winter solstice (supposed to be always 21 December) 105 101 zdsws = 11 + nday_year … … 113 109 DO ji = 1, jpi 114 110 ztmp = zconvrad * gphit(ji,jj) 115 zaaa(ji,jj) = SIN(ztmp) * SIN(zdecrad)116 zbbb(ji,jj) = COS(ztmp) * COS(zdecrad)111 raa(ji,jj) = SIN( ztmp ) * SIN( zdecrad ) 112 rbb(ji,jj) = COS( ztmp ) * COS( zdecrad ) 117 113 END DO 118 114 END DO … … 120 116 ! Compute the time of dawn and dusk 121 117 122 ! zab to test if the day time is equal to 0, less than 24h of full day123 zab(:,:) = -zaaa(:,:) / zbbb(:,:)118 ! rab to test if the day time is equal to 0, less than 24h of full day 119 rab(:,:) = -raa(:,:) / rbb(:,:) 124 120 DO jj = 1, jpj 125 121 DO ji = 1, jpi 126 IF ( ABS(zab(ji,jj)) < 1 ) THEN 127 ! day duration is less than 24h 122 IF ( ABS(rab(ji,jj)) < 1 ) THEN ! day duration is less than 24h 128 123 ! When is it night? 129 ztx = zinvtwopi * (ACOS( zab(ji,jj)) - zccc(ji,jj))130 ztest = - zbbb(ji,jj) * SIN( zccc(ji,jj) + ztwopi * ztx )124 ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 125 ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx ) 131 126 ! is it dawn or dusk? 132 127 IF ( ztest > 0 ) THEN 133 zdawn(ji,jj) = ztx134 zdusk(ji,jj) = ztmd(ji,jj) + ( ztmd(ji,jj) - zdawn(ji,jj) )128 rdawn(ji,jj) = ztx 129 rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) ) 135 130 ELSE 136 zdusk(ji,jj) = ztx137 zdawn(ji,jj) = ztmd(ji,jj) - ( zdusk(ji,jj) - ztmd(ji,jj) )131 rdusk(ji,jj) = ztx 132 rdawn(ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) ) 138 133 ENDIF 139 134 ELSE 140 zdawn(ji,jj) = ztmd(ji,jj) + 0.5141 zdusk(ji,jj) = zdawn(ji,jj)135 rdawn(ji,jj) = rtmd(ji,jj) + 0.5 136 rdusk(ji,jj) = rdawn(ji,jj) 142 137 ENDIF 143 138 END DO 144 139 END DO 145 zdawn(:,:) = MOD((zdawn(:,:) + 1.), 1.)146 zdusk(:,:) = MOD((zdusk(:,:) + 1.), 1.)140 rdawn(:,:) = MOD((rdawn(:,:) + 1.), 1.) 141 rdusk(:,:) = MOD((rdusk(:,:) + 1.), 1.) 147 142 148 143 149 144 ! 2.2 Compute the scalling function: 150 145 ! S* = the inverse of the time integral of the diurnal cycle from dawm to dusk 151 152 146 DO jj = 1, jpj 153 147 DO ji = 1, jpi 154 IF ( ABS(zab(ji,jj)) < 1 ) THEN 155 ! day duration is less than 24h 156 IF ( zdawn(ji,jj) < zdusk(ji,jj) ) THEN 157 ! day time in one part 158 zscal(ji,jj) = fintegral(zdawn(ji,jj), zdusk(ji,jj), zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 159 zscal(ji,jj) = 1. / zscal(ji,jj) 160 ELSE 161 ! day time in two parts 162 zscal(ji,jj) = fintegral(0., zdusk(ji,jj), zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) & 163 & + fintegral(zdawn(ji,jj), 1., zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 164 zscal(ji,jj) = 1. / zscal(ji,jj) 148 IF ( ABS(rab(ji,jj)) < 1 ) THEN ! day duration is less than 24h 149 IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN ! day time in one part 150 rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 151 rscal(ji,jj) = 1. / rscal(ji,jj) 152 ELSE ! day time in two parts 153 rscal(ji,jj) = fintegral(0., rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) & 154 & + fintegral(rdawn(ji,jj), 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 155 rscal(ji,jj) = 1. / rscal(ji,jj) 165 156 ENDIF 166 157 ELSE 167 IF ( zaaa(ji,jj) > zbbb(ji,jj) ) THEN 168 ! 24h day 169 zscal(ji,jj) = fintegral(0., 1., zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 170 zscal(ji,jj) = 1. / zscal(ji,jj) 171 ELSE 172 ! No day 173 zscal(ji,jj) = 0. 158 IF ( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 159 rscal(ji,jj) = fintegral(0., 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 160 rscal(ji,jj) = 1. / rscal(ji,jj) 161 ELSE ! No day 162 rscal(ji,jj) = 0.e0 174 163 ENDIF 175 164 ENDIF … … 178 167 ! 179 168 ztmp = rday / rdt 180 zscal(:,:) = zscal(:,:) * ztmp169 rscal(:,:) = rscal(:,:) * ztmp 181 170 182 171 ENDIF 183 172 184 ! 3. compute qsr with the diurnal cycle185 ! ----------------------- 173 ! 3. update qsr with the diurnal cycle 174 ! ------------------------------------ 186 175 187 176 DO jj = 1, jpj 188 177 DO ji = 1, jpi 189 IF ( ABS(zab(ji,jj)) < 1 ) THEN 190 ! day duration is less than 24h 191 IF ( zdawn(ji,jj) < zdusk(ji,jj) ) THEN 192 ! day time in one part 193 zlousd = MAX(zlo, zdawn(ji,jj)) 194 zlousd = MIN(zlousd, zup) 195 zupusd = MIN(zup, zdusk(ji,jj)) 196 zupusd = MAX(zupusd, zlo) 197 ztmp = fintegral(zlousd, zupusd, zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 198 pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * zscal(ji,jj) 199 ELSE 200 ! day time in two parts 201 zlousd = MIN(zlo, zdusk(ji,jj)) 202 zupusd = MIN(zup, zdusk(ji,jj)) 203 ztmp1 = fintegral(zlousd, zupusd, zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 204 zlousd = MAX(zlo, zdawn(ji,jj)) 205 zupusd = MAX(zup, zdawn(ji,jj)) 206 ztmp2 = fintegral(zlousd, zupusd, zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 207 ztmp = ztmp1 + ztmp2 208 pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * zscal(ji,jj) 209 ENDIF 210 ELSE 211 IF ( zaaa(ji,jj) > zbbb(ji,jj) ) THEN 212 ! 24h day 213 ztmp = fintegral(zlo, zup, zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 214 pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * zscal(ji,jj) 215 ELSE 216 ! No day 217 pqsr(ji,jj) = 0. 218 ENDIF 178 IF( ABS(rab(ji,jj)) < 1 ) THEN ! day duration is less than 24h 179 ! 180 IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN ! day time in one part 181 zlousd = MAX(zlo, rdawn(ji,jj)) 182 zlousd = MIN(zlousd, zup) 183 zupusd = MIN(zup, rdusk(ji,jj)) 184 zupusd = MAX(zupusd, zlo) 185 ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 186 pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * rscal(ji,jj) 187 ! 188 ELSE ! day time in two parts 189 zlousd = MIN(zlo, rdusk(ji,jj)) 190 zupusd = MIN(zup, rdusk(ji,jj)) 191 ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 192 zlousd = MAX(zlo, rdawn(ji,jj)) 193 zupusd = MAX(zup, rdawn(ji,jj)) 194 ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 195 ztmp = ztmp1 + ztmp2 196 pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * rscal(ji,jj) 197 ENDIF 198 ELSE ! 24h light or 24h night 199 ! 200 IF( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 201 ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 202 pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * rscal(ji,jj) 203 ! 204 ELSE ! No day 205 pqsr(ji,jj) = 0.e0 206 ENDIF 219 207 ENDIF 220 208 END DO 221 209 END DO 222 210 ! 223 211 END SUBROUTINE sbc_dcy 224 212
Note: See TracChangeset
for help on using the changeset viewer.
