Changeset 55 for trunk/src/SIMULS_IRCAAM/progfiltrage_simulation.f
- Timestamp:
- 02/10/09 15:03:49 (15 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/SIMULS_IRCAAM/progfiltrage_simulation.f
r22 r55 1 2 c PROGRAM Filtrage 3 4 c filtrage sur olr JAN-DEC de 1968 5 6 parameter (nb=3538.,period1=10.,period2=30.) 7 8 dimension olr(24,15,nb),vb(nb) 9 dimension vvb(nb),vvvb(nb),olrf(24,15,nb) 10 11 open(1,file='olr_ctl.dat' 12 *,form='unformatted',access='direct',recl=nb*24*15*4) 13 read(1,rec=1) 14 *(((olr(i,j,k),i=1,24),j=1,15),k=1,nb) 15 close(1) 1 PROGRAM progfiltrage_simulation 2 C+ 3 C 4 C NAME 5 C ==== 6 C 7 C ``progfiltrage_simulation.x`` 8 C 9 C SYNOPSIS 10 C ======== 11 C 12 C ``progfiltrage_simulation.x`` 13 C 14 C DESCRIPTION 15 C =========== 16 C 17 C filtrage sur olr JAN-DEC de 1968 18 C ++ 19 C 20 C From olr_nofiltre_arpege_ \ *simulation*\ .dat. 21 C ``progfiltrage_simulation.x`` compute ++. 22 C 23 C ``progfiltrage_simulation.x`` write 24 C olr_filtre_\ *period1*\ -\ *period2*\ d_arpege_ \ *simulation*\ .dat. 25 C 26 C 27 C CAUTIONS 28 C ======== 29 C 30 C On MAC ppc, this program must be compile with ``-fendian=little`` 31 C for g95 32 C according to machineformat used in prepare_olr_filtre_simulation.m_. 33 C 34 C EXAMPLES 35 C ======== 36 C 37 C Following line read ``olr_nofiltre_arpege_AfNQIVIV.dat``, 38 C compute ++, 39 C write ``olr_filtre_10-30d_arpege_AfNQIVIV.dat`` 40 C :: 41 C 42 C $ echo "arpege AfNQIVIV 10 30" | progfiltrage_simulation.x - 43 C 44 C SEE ALSO 45 C ======== 46 C 47 C prepare_olr_filtre_simulation.m_ 48 C 49 C .. _prepare_olr_filtre_simulation.m : prepare_olr_filtre_simulation.m.html 50 C 51 C forfilter.f_ 52 C 53 C .. _forfilter.f : forfilter.f.html 54 C 55 C TODO 56 C ==== 57 C 58 C underflow tracking (see ifort failure) 59 C 60 C REAL ou DOUBLE PRECISION dans les fichiers matlab 61 C 62 C improve description 63 C 64 C improve file pb (err in open/read/write) 65 C 66 C F90 modules 67 C 68 C portability check because of getenv : now (20090203) written for 69 C working on zeus using g95 (GCC 4.0.3 (g95!) Jun 18 2006) 70 C 71 C use PT fortran tool box ?! when public ... 72 C 73 C direct access file ? seems to work but see also 74 C http://objectmix.com/fortran/305264-read-binaries-matlab.html 75 C 76 C EVOLUTIONS 77 C ========== 78 C 79 C $Id$ 80 C 81 C - fplod 2009-02-10T13:57:29Z aedon.locean-ipsl.upmc.fr (Darwin) 82 C 83 C * bug fix on output file name to be compatible to convention 84 C (not fixed period number of digits) 85 C 86 C - fplod 2009-02-02T13:45:57Z aedon.locean-ipsl.upmc.fr (Darwin) 87 C 88 C * created from progfiltrage_10_30.f to replace it as well as 89 C progfiltrage_30_100.f 90 C 91 C * Subroutines externalized in forfilter.f 92 C 93 C * add parameters ircaam_dataset and simulation according to 94 C prepare_olr_filtre_simulation.m 95 C 96 C * consolidation : implicit none, inquire iolength, force 97 C LITTLE ENDIAN 98 C 99 C- 100 C 101 IMPLICIT NONE 102 103 EXTERNAL FILTRE 104 105 INTEGER, PARAMETER :: nb = 3660 106 INTEGER :: period1 107 INTEGER :: period2 108 C 109 INTEGER :: i 110 INTEGER :: j 111 INTEGER :: k 112 INTEGER :: l 113 INTEGER :: lon0 114 INTEGER :: lat0 115 C 116 REAL :: olr(24,15,nb),vb(nb) 117 REAL :: vvb(nb),vvvb(nb),olrf(24,15,nb) 118 C 119 CHARACTER (LEN=255) :: ircaam_id 120 CHARACTER (LEN=255) :: ircaam_od 121 C 122 CHARACTER (LEN=255) :: fullfilename 123 CHARACTER (LEN=255) :: fullfilename_format 124 INTEGER :: record_length = 0 125 C 126 CHARACTER (LEN=6) :: ircaam_dataset 127 CHARACTER (LEN=8) :: simulation 128 LOGICAL :: vsel = .TRUE. 129 130 write (*,*) 'enter ircaam_dataset simulation period1 period2' 131 read (*,*) ircaam_dataset,simulation,period1,period2 132 133 write (*,*) 'ircaam_dataset : ', ircaam_dataset 134 write (*,*) 'simulation : ', simulation 135 write (*,*) 'period1 : ', period1 136 write (*,*) 'period2 : ', period2 137 138 C validation of ircaam_dataset 139 IF (ircaam_dataset /= 'arpege') THEN 140 vsel = .FALSE. 141 WRITE(*,*)' eee : pb validation ircaam_dataset' 142 ENDIF 143 144 C validation of simulation 145 IF ((simulation /= 'AfNQIVIV') .AND. 146 + (simulation /= 'TrNQIVIV') .AND. 147 + (simulation /= 'AsNQIVIV') .AND. 148 + (simulation /= 'CtIV') .AND. 149 + (simulation /= 'CtCl ')) THEN 150 vsel = .FALSE. 151 WRITE(*,*)' eee : pb validation simulation' 152 ENDIF 153 154 C validation of period1 and period2 155 156 IF (period1 >= period2) THEN 157 vsel = .FALSE. 158 WRITE(*,*)' eee : pb validation period1 vs period2' 159 ENDIF 160 C 161 IF (.NOT. vsel) THEN 162 GOTO 994 163 ENDIF 164 165 CALL GETENV('IRCAAM_ID',ircaam_id) 166 C++ ircaam_id='/usr/temp/fplod/ircaam_d/' 167 CALL GETENV('IRCAAM_OD',ircaam_od) 168 C++ ircaam_od='/usr/temp/fplod/ircaam_d/' 169 170 WRITE(fullfilename,'(A,A,A,A,A,A)') 171 +TRIM(ircaam_id), 172 +'olr_nofiltre_', 173 +TRIM(ircaam_dataset), 174 +'_', 175 +TRIM(simulation), 176 +'.dat' 177 178 C ask for record length 179 INQUIRE (IOLENGTH=record_length) 180 +(((olr(i,j,k),i=1,24),j=1,15),k=1,nb) 181 182 open(1,file=fullfilename 183 +,status='old' 184 +,form='unformatted',access='direct',recl=record_length) 185 186 read(1,rec=1) 187 +(((olr(i,j,k),i=1,24),j=1,15),k=1,nb) 188 189 CLOSE(1) 16 190 17 191 c l'ordre maximal du filtrage est tel que 2*MOR+1 < nb … … 19 193 20 194 do lon0=1,24 21 do lat0=1,15 22 do i=1,nb 23 vb(i)=olr(lon0,lat0,i) 195 do lat0=1,15 196 do i=1,nb 197 vb(i)=olr(lon0,lat0,i) 198 enddo 199 call filtre(vb,vvb,REAL(period1),nb) 200 call filtre(vb,vvvb,REAL(period2),nb) 201 do i=1,nb 202 olrf(lon0,lat0,i)=vvb(i)-vvvb(i) 203 enddo 204 enddo 24 205 enddo 25 call filtre(vb,vvb,period1,nb) 26 call filtre(vb,vvvb,period2,nb) 27 do i=1,nb 28 olrf(lon0,lat0,i)=vvb(i)-vvvb(i) 29 enddo 30 enddo 31 enddo 32 c print*,(olrf(1,1,j),j=1,nb) 33 34 open(2,file='olrf10-30_clim_sst.80.dat' 35 *,form='unformatted',access='direct',recl=24*15*nb*4) 36 write(2,rec=1)(((olrf(lon0,lat0,l),lon0=1,24), 37 * lat0=1,15),l=1,nb) 206 C print*,(olrf(1,1,j),j=1,nb) 207 208 C adjust format of output filename according to period1 and period2 values 209 fullfilename_format = '(A,A,' 210 IF (period1 < 10 ) THEN 211 fullfilename_format = TRIM(fullfilename_format)//'I1,' 212 ENDIF 213 IF ((period1 >= 10) .AND. (period1 < 100)) THEN 214 fullfilename_format = TRIM(fullfilename_format)//'I2,' 215 ENDIF 216 IF (period1 >= 100 ) THEN 217 fullfilename_format = TRIM(fullfilename_format)//'I3,' 218 ENDIF 219 fullfilename_format = TRIM(fullfilename_format) //'A,' 220 IF (period2 < 10 ) THEN 221 fullfilename_format = TRIM(fullfilename_format)//'I1,' 222 ENDIF 223 IF ((period2 >= 10) .AND. (period2 < 100)) THEN 224 fullfilename_format = TRIM(fullfilename_format)//'I2,' 225 ENDIF 226 IF (period2 >= 100 ) THEN 227 fullfilename_format = TRIM(fullfilename_format)//'I3,' 228 ENDIF 229 230 fullfilename_format = TRIM(fullfilename_format)//'A,A,A,A,A)' 231 232 WRITE(fullfilename,fullfilename_format) 233 +TRIM(ircaam_od), 234 +'olr_filtre_', 235 +period1, 236 +'-', 237 +period2, 238 +'d_', 239 +TRIM(ircaam_dataset), 240 +'_', 241 +TRIM(simulation), 242 +'.dat' 243 244 open(2,file=fullfilename 245 +,status='replace' 246 +,form='unformatted',access='direct',recl=24*15*nb*4) 247 248 write(2,rec=1)(((olrf(lon0,lat0,l),lon0=1,24), 249 + lat0=1,15),l=1,nb) 38 250 close(2) 39 251 40 END 41 42 SUBROUTINE FILTRE(F,F1,PERIO,N) 43 PARAMETER(KOR=4,JOR=4,MOR=50) 44 DIMENSION F(N),W(-MOR:mor),G(-MOR:mor),F1(N) 45 PI=ACOS(-1.) 46 FC=1./PERIO 47 CALL KISER(G,MOR) 48 DO 1 I=-MOR,MOR 49 IF (I.EQ.0) THEN 50 W(I)=2.*FC 51 ELSE 52 W(I)=SIN(2.*PI*FC*FLOAT(I))/(PI*FLOAT(I))*G(I) 53 ENDIF 54 1 CONTINUE 55 56 DO 2 I=1,N 57 F1(I)=0. 58 AT=0. 59 60 IF(I.LE.KOR) THEN 61 L1=-KOR 62 L2=I-1 63 ENDIF 64 65 IF((I.GE.KOR+1).AND.(I.LE.MOR)) THEN 66 L1=-I+1 67 L2=I-1 68 ENDIF 69 70 IF((I.GE.MOR+1).AND.(I.LE.N-MOR)) THEN 71 L1=-MOR 72 L2=MOR 73 ENDIF 74 75 IF((I.GE.N-MOR+1).AND.(I.LE.N-JOR)) THEN 76 L1=-N+I 77 L2=N-I 78 ENDIF 79 80 IF(I.GE.N-JOR+1) THEN 81 L1=-N+I 82 L2=JOR 83 ENDIF 84 85 DO 3 K=L1,L2 86 F1(I)=F1(I)+W(K)*F(I-K) 87 AT=AT+W(K) 88 3 CONTINUE 89 F1(I)=F1(I)/AT 90 2 CONTINUE 91 92 RETURN 252 994 CONTINUE 253 STOP 93 254 END 94 95 SUBROUTINE KISER(W,MOR)96 PARAMETER(LOR=100)97 DIMENSION W(-MOR:MOR),CO(-LOR:LOR)98 99 A=30.100 101 IF(A.LE.21.) THEN102 ALPHA=0.103 ENDIF104 105 IF((A.LT.50.).AND.(A.GT.21.)) THEN106 ALPHA=0.5842*(A-21.)**0.4+0.07886*(A-21.)107 ENDIF108 109 IF(A.GE.50.) THEN110 ALPHA=0.1102*(A-8.7)111 ENDIF112 113 DO 2 I=-MOR,MOR114 CO(I)=ALPHA*SQRT(1.-(FLOAT(I)/FLOAT(MOR))**2)115 116 N=0117 AS=1.118 AU=1.119 AS1=1.120 AU1=1.121 122 DO 10 K=1,200123 N=N+1124 125 AU=AU*((CO(I)/2.)/FLOAT(N))**2126 AS=AS+AU127 128 AU1=AU1*((ALPHA/2.)/FLOAT(N))**2129 AS1=AS1+AU1130 10 CONTINUE131 W(I)=AS/AS1132 2 CONTINUE133 RETURN134 END135 136
Note: See TracChangeset
for help on using the changeset viewer.