Changeset 6370
- Timestamp:
- 2016-03-08T10:20:12+01:00 (8 years ago)
- Location:
- branches/UKMO/dev_5518_tide_analysis_restart/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 4 edited
-
BDY/bdytides.F90 (modified) (3 diffs)
-
DIA/diaharmana.F90 (modified) (13 diffs)
-
DOM/istate.F90 (modified) (1 diff)
-
step.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_5518_tide_analysis_restart/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
r6116 r6370 52 52 INTEGER, PUBLIC, PARAMETER :: jptides_max = 15 !: Max number of tidal contituents 53 53 LOGICAL, PUBLIC :: ln_harm_ana !: =T Compute harmonic Analysis 54 LOGICAL, PUBLIC :: ln_harmana_read !: =T Decide to do the analysis 55 !from scratch or continue previous run 54 56 55 57 !$AGRIF_DO_NOT_TREAT … … 92 94 TYPE(MAP_POINTER), DIMENSION(jpbgrd) :: ibmap_ptr !: array of pointers to nbmap 93 95 !! 94 NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj, ln_harm_ana 96 NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj, ln_harm_ana, ln_harmana_read 95 97 !!---------------------------------------------------------------------- 96 98 … … 128 130 IF(lwp) WRITE(numout,*) ' Number of tidal components to read: ', nb_harmo 129 131 IF(lwp) WRITE(numout,*) ' Use PCOMS harmonic ananalysis or not: ', ln_harm_ana 132 IF(lwp) WRITE(numout,*) ' Read in previous days harmonic data or start afresh: ', ln_harmana_read 130 133 IF(lwp) THEN 131 134 WRITE(numout,*) ' Tidal components: ' -
branches/UKMO/dev_5518_tide_analysis_restart/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90
r6117 r6370 25 25 USE daymod ! calendar 26 26 USE tideini 27 USE restart 28 USE ioipsl, ONLY : ju2ymds ! for calendar 29 27 30 IMPLICIT NONE 28 31 PRIVATE … … 41 44 REAL(wp), DIMENSION(nharm_max) :: & 42 45 om_tide ! tidal frequencies ( rads/sec) 43 REAL(wp), DIMENSION(nhm_max) :: &44 bzz (nhm_max),c(nhm_max),x(nhm_max)! work arrays46 REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:) :: & 47 bzz,c,x ! work arrays 45 48 REAL(wp) :: cca,ssa,zm,bt 46 49 REAL(wp) :: ccau,ccav,ssau,ssav 47 REAL(wp), PUBLIC, DIMENSION(15) :: anau, anav, anaf ! nodel/phase corrections used by diaharmana 48 ! 49 REAL(wp),ALLOCATABLE,SAVE, DIMENSION(:,:,:,:,:) :: & 50 bz ! work arrays for u,v anaylsis 51 REAL(wp),ALLOCATABLE,SAVE, DIMENSION(:,:,:,:,:) :: & 52 cosamp, & 53 sinamp 50 REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf ! nodel/phase corrections used by diaharmana 54 51 ! 55 52 REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: & 56 bssh ! work array for ssh anaylsis53 bssh,bubar,bvbar ! work array for ssh anaylsis 57 54 REAL(wp), ALLOCATABLE,SAVE, DIMENSION(:,:,:) :: & 58 cosampz, &59 sinampz 55 cosampz, cosampu,cosampv, & 56 sinampz, sinampu,sinampv 60 57 REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a 61 58 ! … … 64 61 REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) :: & 65 62 huout,hvout,guout,gvout ! arrays for output 63 REAL(wp), PUBLIC :: fjulday_startharm !: Julian Day since start of harmonic analysis 66 64 67 65 !! * Substitutions … … 105 103 !! * local declarations 106 104 INTEGER :: ji, jj, jk, & ! dummy loop arguments 107 ih,i1,i2,iv 105 ih,i1,i2,iv,jgrid 106 108 107 109 108 !!-------------------------------------------------------------------- … … 116 115 !nharm=ncpt_anal 117 116 nharm=nb_harmo 118 nhm=2*n harm+1117 nhm=2*nb_harmo+1 119 118 120 119 c(1)=1.0 121 120 122 do ih=1,n harm123 c(2*ih)= cos(adatrj*86400._wp*om_tide(ih))124 c(2*ih+1)= sin(adatrj*86400._wp*om_tide(ih))121 do ih=1,nb_harmo 122 c(2*ih)= cos((fjulday-fjulday_startharm)*86400._wp*om_tide(ih) ) 123 c(2*ih+1)= sin((fjulday-fjulday_startharm)*86400._wp*om_tide(ih) ) 125 124 enddo 126 125 … … 129 128 do jj=1,jpj 130 129 do ih=1,nhm 131 bssh(ih,ji,jj)=bssh(ih,ji,jj)+c(ih)*sshn(ji,jj) 130 bssh (ih,ji,jj)=bssh (ih,ji,jj)+c(ih)*sshn(ji,jj) 131 bubar(ih,ji,jj)=bubar(ih,ji,jj)+c(ih)*un_b(ji,jj) 132 bvbar(ih,ji,jj)=bvbar(ih,ji,jj)+c(ih)*vn_b(ji,jj) 132 133 enddo 133 134 enddo 134 135 enddo 135 136 ! 136 do ji=1,jpi 137 do jj=1,jpj 138 do jk=1,jpk 139 do ih=1,(2*nharm+1) 140 bz(1,ih,ji,jj,jk)=bz(1,ih,ji,jj,jk)+c(ih)*un(ji,jj,jk) !ua(ji,jj,jk) 141 bz(2,ih,ji,jj,jk)=bz(2,ih,ji,jj,jk)+c(ih)*vn(ji,jj,jk) !va(ji,jj,jk) 142 enddo 143 enddo 144 enddo 145 enddo 146 IF(lwp) WRITE(666,'(4E15.5)') adatrj*86400._wp, sshn(10,10),bssh(2,10,10),bssh(3,10,10) 137 ! IF(lwp) WRITE(666,'(4E15.5)') adatrj*86400._wp, sshn(10,10),bssh(2,10,10),bssh(3,10,10) 147 138 ! 148 139 do i1=1,nhm … … 159 150 IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run' 160 151 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 161 162 ! 163 cosamp=0.0_wp 164 sinamp=0.0_wp 152 CALL harm_rst_write(kt) ! Dump out data for a restarted run 153 154 155 156 ! 157 cosampu=0.0_wp 158 sinampu=0.0_wp 159 cosampv=0.0_wp 160 sinampv=0.0_wp 165 161 cosampz=0.0_wp 166 162 sinampz=0.0_wp 167 !velocities 168 do ji=1,jpi 169 do jj=1,jpj 170 do iv=1,nvab 171 do jk=1,jpk 172 bt=1.0 173 do ih=1,nhm 174 bzz(ih)=bz(iv,ih,ji,jj,jk) 175 bt=bt*bzz(ih) 176 enddo 177 ! 178 do i1=1,nhm 179 do i2=1,nhm 180 a(i1,i2)=cc(i1,i2) 181 enddo 182 enddo 183 ! 184 ! now do gaussian elimination of the system 185 ! a * x = b 186 ! the matrix x is (a0,a1,b1,a2,b2 ...) 187 ! the matrix a and rhs b solved here for x 188 189 if(bt.ne.0.) then 190 call gelim(a,bzz,x,nhm) 191 ! 192 do ih=1,nharm 193 cosamp(iv,ih,ji,jj,jk)=x(ih*2) 194 sinamp(iv,ih,ji,jj,jk)=x(ih*2+1) 195 enddo 196 cosamp(iv,0,ji,jj,jk)=x(1) 197 ! 198 endif 199 ! 200 enddo 201 ! next variable 202 enddo 203 enddo 204 enddo 205 ! 206 ! surface elevation 163 ! 164 do jgrid=1,3 ! elevation, Ubar,Vbar 207 165 do ji=1,jpi 208 166 do jj=1,jpj 209 167 bt=1.0 210 168 do ih=1,nhm 211 bzz(ih)=bssh(ih,ji,jj) 169 if(jgrid .eq. 1) then 170 bzz(ih)=bssh(ih,ji,jj) 171 endif 172 if(jgrid .eq. 2) then 173 bzz(ih)=bubar(ih,ji,jj) 174 endif 175 if(jgrid .eq. 3) then 176 bzz(ih)=bvbar(ih,ji,jj) 177 endif 212 178 bt=bt*bzz(ih) 213 179 enddo … … 228 194 call gelim(a,bzz,x,nhm) 229 195 230 do ih=1,nharm 231 cosampz(ih,ji,jj)=x(ih*2) 232 sinampz(ih,ji,jj)=x(ih*2+1) 196 do ih=1,nb_harmo 197 if(jgrid .eq. 1) then 198 cosampz(ih,ji,jj)=x(ih*2) 199 sinampz(ih,ji,jj)=x(ih*2+1) 200 endif 201 if(jgrid .eq. 2) then 202 cosampu(ih,ji,jj)=x(ih*2) 203 sinampu(ih,ji,jj)=x(ih*2+1) 204 endif 205 if(jgrid .eq. 3) then 206 cosampv(ih,ji,jj)=x(ih*2) 207 sinampv(ih,ji,jj)=x(ih*2+1) 208 endif 233 209 enddo 210 if(jgrid .eq. 1) then 234 211 cosampz(0,ji,jj)=x(1) 235 212 sinampz(0,ji,jj)=0.0_wp 213 endif 214 if(jgrid .eq. 2) then 215 cosampu(0,ji,jj)=x(1) 216 sinampu(0,ji,jj)=0.0_wp 217 endif 218 if(jgrid .eq. 3) then 219 cosampv(0,ji,jj)=x(1) 220 sinampv(0,ji,jj)=0.0_wp 221 endif 236 222 endif 237 223 enddo 238 224 enddo 239 if(lwp) write(numout,*) 'cos',cosampz(1,10,10),sinampz(1,10,10) 240 if(lwp) then 241 do i1=1,nhm 242 do i2=1,nhm 243 ! write(numout,*) i1,i2,cc(i1,i2) 244 enddo 245 enddo 246 endif 225 enddo ! jgrid 226 ! 247 227 ! 248 228 CALL harm_ana_out ! output analysis (last time step) 229 249 230 ENDIF 250 231 END SUBROUTINE harm_ana … … 273 254 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 274 255 ! DO ALLOCATIONS 275 ALLOCATE( bz(nvab,nhm_max,jpi,jpj,jpk)) 276 ALLOCATE( cosamp(nvab,0:nharm_max,jpi,jpj,jpk)) 277 ALLOCATE( sinamp(nvab,0:nharm_max,jpi,jpj,jpk) ) 278 ALLOCATE( bssh(nhm_max,jpi,jpj) ) 279 ALLOCATE( cosampz(0:nharm_max,jpi,jpj)) 280 ALLOCATE( sinampz(0:nharm_max,jpi,jpj)) 281 ALLOCATE( cc(nhm_max,nhm_max) ) 282 ALLOCATE( a(nhm_max,nhm_max) ) 256 257 ALLOCATE( bubar(nb_harmo*2+1,jpi,jpj) ) 258 ALLOCATE( bvbar(nb_harmo*2+1,jpi,jpj) ) 259 ALLOCATE( bssh (nb_harmo*2+1,jpi,jpj) ) 260 261 ALLOCATE( cosampz(0:nb_harmo*2+1,jpi,jpj)) 262 ALLOCATE( cosampu(0:nb_harmo*2+1,jpi,jpj)) 263 ALLOCATE( cosampv(0:nb_harmo*2+1,jpi,jpj)) 264 265 ALLOCATE( sinampz(0:nb_harmo*2+1,jpi,jpj)) 266 ALLOCATE( sinampu(0:nb_harmo*2+1,jpi,jpj)) 267 ALLOCATE( sinampv(0:nb_harmo*2+1,jpi,jpj)) 268 269 ALLOCATE( cc(nb_harmo*2+1,nb_harmo*2+1) ) 270 271 ALLOCATE( a(nb_harmo*2+1,nb_harmo*2+1) ) 272 273 ALLOCATE( anau(nb_harmo) ) 274 ALLOCATE( anav(nb_harmo) ) 275 ALLOCATE( anaf(nb_harmo) ) 276 277 ALLOCATE( bzz(nb_harmo*2+1) ) 278 ALLOCATE( x(nb_harmo*2+1) ) 279 ALLOCATE( c(nb_harmo*2+1) ) 280 283 281 ALLOCATE( gout(jpi,jpj)) 284 282 ALLOCATE( hout(jpi,jpj)) 283 285 284 ALLOCATE( guout(jpi,jpj)) 286 285 ALLOCATE( huout(jpi,jpj)) 286 287 287 ALLOCATE( gvout(jpi,jpj)) 288 288 ALLOCATE( hvout(jpi,jpj)) 289 289 290 ! END ALLOCATE 290 291 291 ! do ih=1,ncpt_anal292 292 do ih=1,nb_harmo 293 ! om_tide(ih)= omega_tide(ih) *rpi/(3600._wp*180._wp)294 293 om_tide(ih)= omega_tide(ih) 295 ! WRITE(numout,*) " ih, om_tide(ih):",ih," ",om_tide(ih)," omega_tide(ih):",omega_tide(ih)296 294 enddo 297 bz(:,:,:,:,:)=0.0_wp 298 bssh(:,:,:)=0.0_wp 299 cc=0.0_wp 300 301 anau(:) = 0. 302 anav(:) = 0. 303 anaf(:) = 0. 295 if(ln_harmana_read) then 296 WRITE(numout,*) "Reading previous harmonic data from previous run" 297 ! Need to read in bssh bz, cc anau anav and anaf 298 call harm_rst_read ! This reads in from the previous day 299 ! Currrently the data in in assci format 300 else 301 WRITE(numout,*) "Starting harmonic analysis from Fresh " 302 bssh(:,:,:) = 0.0_wp 303 bubar(:,:,:) = 0.0_wp 304 bvbar(:,:,:) = 0.0_wp 305 306 cc=0.0_wp 307 308 anau(:) = 0.0_wp 309 anav(:) = 0.0_wp 310 anaf(:) = 0.0_wp 311 bzz(:) = 0.0_wp 312 x(:) = 0.0_wp 313 c(:) = 0.0_wp 314 304 315 DO ih = 1, nb_harmo 305 316 anau(ih) = utide(ih) 306 317 anav(ih) = v0tide(ih) 307 318 anaf(ih) = ftide(ih) 308 ! WRITE(numout,*) " EOD anaf(",ih,"):",anaf(ih) ," anau(",ih,"):",anau(ih) , " anav(",ih,"):",anav(ih)309 319 310 320 END DO 321 fjulday_startharm=fjulday !Set this at very start and store 322 endif 311 323 312 324 END SUBROUTINE harm_ana_init … … 332 344 ! 333 345 integer :: n 334 REAL(WP) :: b(n hm_max),a(nhm_max,nhm_max)335 REAL(WP) :: x(n hm_max)346 REAL(WP) :: b(nb_harmo*2+1),a(nb_harmo*2+1,nb_harmo*2+1) 347 REAL(WP) :: x(nb_harmo*2+1) 336 348 INTEGER :: row,col,prow,pivrow,rrow,ntemp 337 349 REAL(WP) :: atemp 338 350 REAL(WP) :: pivot 339 351 REAL(WP) :: m 340 REAL(WP) :: tempsol(n hm_max)352 REAL(WP) :: tempsol(nb_harmo*2+1) 341 353 342 354 … … 406 418 407 419 !! * local declarations 408 INTEGER :: ji, jj, jk, jit, ih, jjk ! dummy loop indices 409 ! INTEGER :: & 410 ! iimi, iima, ipk, it, & ! temporary integers 411 ! ijmi, ijma ! " " 420 INTEGER :: ji, jj, jk, jgrid,jit, ih, jjk ! dummy loop indices 412 421 INTEGER :: nh_T 413 422 INTEGER :: nid_harm 414 423 CHARACTER (len=40) :: & 415 clhstnam , clop1, clop2 ! temporary names424 clhstnamt, clop1, clop2 ! temporary names 416 425 CHARACTER (len=40) :: & 417 clhstunamt, clhstunamm, clhstunamb ! temporary names 418 CHARACTER (len=40) :: & 419 clhstvnamt, clhstvnamm, clhstvnamb ! temporary names 426 clhstnamu, clhstnamv ! temporary names 420 427 REAL(wp) :: & 421 428 zsto1, zsto2, zout, zmax, & ! temporary scalars 422 429 zjulian, zdt, zmdi 423 430 424 !!---------------------------------------------------------------------- 425 ! Define indices of the horizontal output zoom and vertical limit storage 426 ! iimi = 1 ; iima = jpi 427 ! ijmi = 1 ; ijma = jpj 428 ! ipk = jpk 429 ! 430 clhstnam=TRIM(cexper)//'.harm_ana.grid_T' 431 clhstunamt=TRIM(cexper)//'.harm_ana.grid_U_TOP' 432 clhstunamm=TRIM(cexper)//'.harm_ana.grid_U_MID' 433 clhstunamb=TRIM(cexper)//'.harm_ana.grid_U_BOT' 434 clhstvnamt=TRIM(cexper)//'.harm_ana.grid_V_TOP' 435 clhstvnamm=TRIM(cexper)//'.harm_ana.grid_V_MID' 436 clhstvnamb=TRIM(cexper)//'.harm_ana.grid_V_BOT' 437 write(6,*) nproc,nimpp,nimpp+jpi-1,nharm 438 write (clhstnam,'(a,''_'',i3.3)') trim(clhstnam),nproc 439 write (clhstunamt,'(a,''_'',i3.3)') trim(clhstunamt),nproc 440 write (clhstunamm,'(a,''_'',i3.3)') trim(clhstunamm),nproc 441 write (clhstunamb,'(a,''_'',i3.3)') trim(clhstunamb),nproc 442 write (clhstvnamt,'(a,''_'',i3.3)') trim(clhstvnamt),nproc 443 write (clhstvnamm,'(a,''_'',i3.3)') trim(clhstvnamm),nproc 444 write (clhstvnamb,'(a,''_'',i3.3)') trim(clhstvnamb),nproc 445 ! 446 open(66,file=clhstnam) 447 open(67,file=clhstunamt) 448 open(68,file=clhstunamm) 449 open(69,file=clhstunamb) 450 open(70,file=clhstvnamt) 451 open(71,file=clhstvnamm) 452 open(72,file=clhstvnamb) 453 write(6,*) 'out to:',clhstnam 454 455 write(66,*) nimpp,nimpp+jpi-1 456 write(66,*) njmpp,njmpp+jpj-1 457 write(67,*) nimpp,nimpp+jpi-1 458 write(67,*) njmpp,njmpp+jpj-1 459 write(68,*) nimpp,nimpp+jpi-1 460 write(68,*) njmpp,njmpp+jpj-1 461 write(69,*) nimpp,nimpp+jpi-1 462 write(69,*) njmpp,njmpp+jpj-1 463 write(70,*) nimpp,nimpp+jpi-1 464 write(70,*) njmpp,njmpp+jpj-1 465 write(71,*) nimpp,nimpp+jpi-1 466 write(71,*) njmpp,njmpp+jpj-1 467 write(72,*) nimpp,nimpp+jpi-1 468 write(72,*) njmpp,njmpp+jpj-1 469 ! 470 do ih=1,nharm 431 do jgrid=1,3 432 do ih=1,nb_harmo 471 433 hout = 0.0 472 434 gout = 0.0 473 435 do jj=1,nlcj 474 436 do ji=1,nlci 475 cca=cosampz(ih,ji,jj) 476 ssa=sinampz(ih,ji,jj) 437 if(jgrid .eq. 1) then ! SSH 438 cca=cosampz(ih,ji,jj) 439 ssa=sinampz(ih,ji,jj) 440 endif 441 if(jgrid .eq. 2) then ! UBAR 442 cca=cosampu(ih,ji,jj) 443 ssa=sinampu(ih,ji,jj) 444 endif 445 if(jgrid .eq. 3) then ! VBAR 446 cca=cosampv(ih,ji,jj) 447 ssa=sinampv(ih,ji,jj) 448 endif 449 477 450 hout(ji,jj)=sqrt(cca**2+ssa**2) 478 451 … … 502 475 endif 503 476 if (gout(ji,jj).ne.0) then 504 gout(ji,jj)=gout(ji,jj)+anau(ih)+anav(ih) 477 !Convert Rad to degree and take 478 !modulus 479 gout(ji,jj)=gout(ji,jj)+MOD( (anau(ih)+anav(ih))/rad , 360.0) 505 480 if (gout(ji,jj).gt.360.0) then 506 481 gout(ji,jj)=gout(ji,jj)-360.0 507 482 else if (gout(ji,jj).lt.0) then 508 gout(ji,jj)=gout(ji,jj)+ 180.0483 gout(ji,jj)=gout(ji,jj)+360.0 509 484 endif 510 485 endif 511 486 enddo 512 487 enddo 513 write(66,'(50e15.5)') hout 514 write(66,'(50e15.5)') gout 488 !NETCDF OUTPUT 489 if(jgrid==1) then 490 WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'amp' 491 WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'phase' 492 CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp', hout(:,:) ) 493 CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase', gout(:,:) ) 494 endif 495 if(jgrid==2) then 496 CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_u', hout(:,:) ) 497 CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_u', gout(:,:) ) 498 endif 499 if(jgrid==3) then 500 CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_v', hout(:,:) ) 501 CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_v', gout(:,:) ) 502 endif 503 enddo 515 504 enddo 516 close (66)517 518 ! output u and v519 do jjk=1,3520 do ih=1,nharm521 huout = 0.0522 guout = 0.0523 hvout = 0.0524 gvout = 0.0525 do jj=1,nlcj526 do ji=1,nlci527 if (jjk.eq.1) then528 jk=1529 else if (jjk.eq.2) then530 jk=max(1,mbathy(ji,jj)/2)531 else532 jk=max(1,mbathy(ji,jj)-1)533 endif534 535 ccau=cosamp(1,ih,ji,jj,jk)536 ssau=sinamp(1,ih,ji,jj,jk)537 huout(ji,jj)=sqrt(ccau**2+ssau**2)538 ccav=cosamp(2,ih,ji,jj,jk)539 ssav=sinamp(2,ih,ji,jj,jk)540 hvout(ji,jj)=sqrt(ccav**2+ssav**2)541 542 ! U543 if (ccau.ne.0.0) then544 guout(ji,jj)=(180.0/rpi)*atan(ssau/ccau)545 else546 if (ssau.ne.0.0) then547 if (ssau.gt. 0.0) then548 guout(ji,jj)=90.0549 else550 guout(ji,jj)=270.0551 endif552 else553 guout(ji,jj)=0.0554 endif555 endif556 if (guout(ji,jj).lt.0) then557 guout(ji,jj)=guout(ji,jj)+180.0558 endif559 if (ssau.lt.0) then560 guout(ji,jj)=guout(ji,jj)+180.0561 endif562 563 if (huout(ji,jj).ne.0) then564 huout(ji,jj)=huout(ji,jj)/anaf(ih)565 endif566 if (guout(ji,jj).ne.0) then567 guout(ji,jj)=guout(ji,jj)+anau(ih)+anav(ih)568 if (guout(ji,jj).gt.360.0) then569 guout(ji,jj)=guout(ji,jj)-360.0570 else if (guout(ji,jj).lt.0) then571 guout(ji,jj)=guout(ji,jj)+180.0572 endif573 endif574 575 ! V576 if (ccav.ne.0.0) then577 gvout(ji,jj)=(180.0/rpi)*atan(ssav/ccav)578 else579 if (ssav.ne.0.0) then580 if (ssav.gt. 0.0) then581 gvout(ji,jj)=90.0582 else583 gvout(ji,jj)=270.0584 endif585 else586 gvout(ji,jj)=0.0587 endif588 endif589 if (gvout(ji,jj).lt.0) then590 gvout(ji,jj)=gvout(ji,jj)+180.0591 endif592 if (ssav.lt.0) then593 gvout(ji,jj)=gvout(ji,jj)+180.0594 endif595 596 if (hvout(ji,jj).ne.0) then597 hvout(ji,jj)=hvout(ji,jj)/anaf(ih)598 endif599 if (gvout(ji,jj).ne.0) then600 gvout(ji,jj)=gvout(ji,jj)+anau(ih)+anav(ih)601 if (gvout(ji,jj).gt.360.0) then602 gvout(ji,jj)=gvout(ji,jj)-360.0603 else if (gvout(ji,jj).lt.0) then604 gvout(ji,jj)=gvout(ji,jj)+180.0605 endif606 endif607 enddo608 enddo609 write(66+jjk,'(50e15.5)') huout610 write(66+jjk,'(50e15.5)') guout611 write(69+jjk,'(50e15.5)') hvout612 write(69+jjk,'(50e15.5)') gvout613 enddo614 enddo615 close(67)616 close(68)617 close(69)618 close(70)619 close(71)620 close(72)621 505 622 ! CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit623 ! & 1,jpi, 1, jpj, &624 ! & 0, zjulian, (nitend-nit000+1)*zdt, nh_T, nid_harm, domain_id=nidom )625 ! CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit626 ! & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &627 ! & 0, zjulian, zdt, nh_T, nid_harm, domain_id=nidom )628 ! CALL histdef( nid_harm, "harm_ana", "Sea Surface Height tidal analysis" , "m" , & ! ssh629 ! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, "inst(x)", (nitend-nit000+1)*zdt, (nitend-nit000+1)*zdt )!630 !631 ! CALL histwrite( nid_harm, "cosampz", it, zw2d , ndim_hT, ndex_hT ) ! sea surface height632 506 ! 633 507 END SUBROUTINE harm_ana_out 508 509 SUBROUTINE harm_rst_write(kt) 510 !!---------------------------------------------------------------------- 511 !! *** ROUTINE harm_ana_init *** 512 !! 513 !! ** Purpose : To write out cummulated Tidal Harmomnic data to file for 514 !! restarting 515 !! 516 !! ** Method : restart files will be dated by default 517 !! 518 !! ** input : 519 !! 520 !! ** Action : ... 521 !! 522 !! history : 523 !! 0.0 ! 01-16 (Enda O'Dea) Original code 524 !! ASSUMES dated file for rose , can change later to be more generic 525 !!---------------------------------------------------------------------- 526 INTEGER, INTENT(in) :: kt ! ocean time-step 527 INTEGER :: iyear, imonth, iday,ih 528 REAL (wp) :: zsec 529 !! 530 CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character 531 CHARACTER(LEN=50) :: clname ! ocean output restart file name 532 CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file 533 CHARACTER(LEN=250) :: clfinal ! full name 534 535 CALL ju2ymds( fjulday , iyear, imonth, iday, zsec) 536 537 WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 538 ! create the file 539 clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana" 540 clpath = TRIM(cn_ocerst_outdir) 541 IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/' 542 WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname 543 544 545 546 547 !! clhstnam=TRIM(cexper)//'.restart_harm_ana' 548 ! Open a file to write the data into 549 write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc 550 open(66,file=TRIM(clfinal)) 551 write(66,'(1e15.5)') cc 552 !These Three which contain the most data can be moved to a regular 553 !restart file 554 555 CALL iom_rstput( kt, nitrst, numrow, 'Mean_bssh' , bssh(1,:,:) ) 556 CALL iom_rstput( kt, nitrst, numrow, 'Mean_bubar' , bubar(1,:,:) ) 557 CALL iom_rstput( kt, nitrst, numrow, 'Mean_bvbar' , bvbar(1,:,:) ) 558 do ih=1,nb_harmo 559 CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos' , bssh (ih*2 , : , : ) ) 560 CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin' , bssh (ih*2+1, : , : ) ) 561 562 CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos' , bubar(ih*2 , : , : ) ) 563 CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin' , bubar(ih*2+1, : , : ) ) 564 565 CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos' , bvbar(ih*2 , : , : ) ) 566 CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin' , bvbar(ih*2+1, : , : ) ) 567 enddo 568 569 write(66,'(1e15.5)') anau 570 write(66,'(1e15.5)') anav 571 write(66,'(1e15.5)') anaf 572 write(66,'(1e25.20)') fjulday_startharm 573 close(66) 574 575 END SUBROUTINE harm_rst_write 576 577 SUBROUTINE harm_rst_read 578 !!---------------------------------------------------------------------- 579 !! *** ROUTINE harm_ana_init *** 580 !! 581 !! ** Purpose : To read in cummulated Tidal Harmomnic data to file for 582 !! restarting 583 !! 584 !! ** Method : 585 !! 586 !! ** input : 587 !! 588 !! ** Action : ... 589 !! 590 !! history : 591 !! 0.0 ! 01-16 (Enda O'Dea) Original code 592 !! ASSUMES dated file for rose , can change later to be more generic 593 !!---------------------------------------------------------------------- 594 CHARACTER(LEN=50) :: clname ! ocean output restart file name 595 CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file 596 CHARACTER(LEN=250) :: clfinal ! full name 597 INTEGER :: ih 598 599 ! create the file 600 clname = "restart_harm_ana" 601 clpath = TRIM(cn_ocerst_indir) 602 IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/' 603 WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname 604 605 606 607 608 !! clhstnam=TRIM(cexper)//'.restart_harm_ana' 609 ! Open a file to read the data ifrom 610 write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc 611 open(66,file=TRIM(clfinal),status='old') 612 read(66,'(1e15.5)') cc 613 !2D regular fields can be read from normal restart this saves space and handy to 614 !view in netcdf format also. 615 CALL iom_get( numror,jpdom_autoglo, 'Mean_bssh' , bssh(1,:,:) ) 616 CALL iom_get( numror,jpdom_autoglo, 'Mean_bubar' , bubar(1,:,:) ) 617 CALL iom_get( numror,jpdom_autoglo, 'Mean_bvbar' , bvbar(1,:,:) ) 618 619 do ih=1,nb_harmo 620 621 CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos' , bssh (ih*2 , : , : ) ) 622 CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin' , bssh (ih*2+1, : , : ) ) 623 624 CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos' , bubar(ih*2 , : , : ) ) 625 CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin' , bubar(ih*2+1, : , : ) ) 626 627 CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos' , bvbar(ih*2 , : , : ) ) 628 CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin' , bvbar(ih*2+1, : , : ) ) 629 enddo 630 read(66,'(1e15.5)') anau 631 read(66,'(1e15.5)') anav 632 read(66,'(1e15.5)') anaf 633 read(66,'(1e25.20)') fjulday_startharm 634 WRITE(numout,*) 'Checking anaf is correct' 635 WRITE(numout,*) anaf 636 WRITE(numout,*) '---end of anaf checking----' 637 638 close(66) 639 640 END SUBROUTINE harm_rst_read 641 634 642 !!====================================================================== 635 643 #else … … 638 646 !!--------------------------------------------------------------------------------- 639 647 CONTAINS 648 SUBROUTINE harm_rst_write(kt) ! Dummy routine 649 END SUBROUTINE harm_rst_write 650 SUBROUTINE harm_rst_read ! Dummy routine 651 END SUBROUTINE harm_rst_read 640 652 SUBROUTINE harm_ana_out ! Dummy routine 641 653 END SUBROUTINE harm_ana_out -
branches/UKMO/dev_5518_tide_analysis_restart/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r6061 r6370 211 211 ! 212 212 DO jk = 1, jpk 213 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) ) & 214 & + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.) ) * tmask(:,:,jk) 213 tsn(:,:,jk,jp_tem) = 10.0 * tmask(:,:,jk) 215 214 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 216 215 END DO -
branches/UKMO/dev_5518_tide_analysis_restart/NEMOGCM/NEMO/OPA_SRC/step.F90
r6116 r6370 347 347 CALL dia_wri_state( 'output.abort', kstp ) 348 348 ENDIF 349 IF( ln_harm_ana ) CALL harm_ana( kstp ) ! Harmonic analysis of tides 349 350 IF( kstp == nit000 ) THEN 350 351 CALL iom_close( numror ) ! close input ocean restart file … … 352 353 IF( lwm.AND.numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice 353 354 ENDIF 354 IF( lrst_oce 355 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 355 356 356 357 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 365 366 ENDIF 366 367 #endif 367 IF( ln_harm_ana ) CALL harm_ana( kstp ) ! Harmonic analysis of tides368 368 ! 369 369 IF( nn_timing == 1 .AND. kstp == nit000 ) CALL timing_reset
Note: See TracChangeset
for help on using the changeset viewer.
