Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r5836 r6140 11 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 12 12 !! 3.4 ! 2012-10 (A. Weaver and K. Mogensen) Fix for direct initialization 13 !! ! 2014-09 (D. Lea) Local calc_date removed use routine from OBS 14 !! ! 2015-11 (D. Lea) Handle non-zero initial time of day 13 15 !!---------------------------------------------------------------------- 14 16 15 17 !!---------------------------------------------------------------------- 16 18 !! asm_inc_init : Initialize the increment arrays and IAU weights 17 !! calc_date : Compute the calendar date YYYYMMDD on a given step18 19 !! tra_asm_inc : Apply the tracer (T and S) increments 19 20 !! dyn_asm_inc : Apply the dynamic (u and v) increments … … 38 39 #endif 39 40 USE sbc_oce ! Surface boundary condition variables. 41 USE diaobs, ONLY: calc_date ! Compute the calendar date on a given step 40 42 41 43 IMPLICIT NONE … … 43 45 44 46 PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights 45 PUBLIC calc_date !: Compute the calendar date YYYYMMDD on a given step46 47 PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments 47 48 PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments … … 87 88 88 89 !! * Substitutions 89 # include "domzgr_substitute.h90"90 90 # include "vectopt_loop_substitute.h90" 91 91 !!---------------------------------------------------------------------- 92 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)92 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 93 93 !! $Id$ 94 94 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 111 111 INTEGER :: iiauper ! Number of time steps in the IAU period 112 112 INTEGER :: icycper ! Number of time steps in the cycle 113 INTEGER :: iitend_date ! Date YYYYMMDDof final time step114 INTEGER :: iitbkg_date ! Date YYYYMMDDof background time step for Jb term115 INTEGER :: iitdin_date ! Date YYYYMMDDof background time step for DI116 INTEGER :: iitiaustr_date ! Date YYYYMMDDof IAU interval start time step117 INTEGER :: iitiaufin_date ! Date YYYYMMDDof IAU interval final time step118 ! 113 REAL(KIND=dp) :: ditend_date ! Date YYYYMMDD.HHMMSS of final time step 114 REAL(KIND=dp) :: ditbkg_date ! Date YYYYMMDD.HHMMSS of background time step for Jb term 115 REAL(KIND=dp) :: ditdin_date ! Date YYYYMMDD.HHMMSS of background time step for DI 116 REAL(KIND=dp) :: ditiaustr_date ! Date YYYYMMDD.HHMMSS of IAU interval start time step 117 REAL(KIND=dp) :: ditiaufin_date ! Date YYYYMMDD.HHMMSS of IAU interval final time step 118 119 119 REAL(wp) :: znorm ! Normalization factor for IAU weights 120 120 REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights (should be equal to one) … … 178 178 icycper = nitend - nit000 + 1 ! Cycle interval length 179 179 180 CALL calc_date( nit000, nitend , ndate0, iitend_date ) ! Date of final time step 181 CALL calc_date( nit000, nitbkg_r , ndate0, iitbkg_date ) ! Background time for Jb referenced to ndate0 182 CALL calc_date( nit000, nitdin_r , ndate0, iitdin_date ) ! Background time for DI referenced to ndate0 183 CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date ) ! IAU start time referenced to ndate0 184 CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date ) ! IAU end time referenced to ndate0 185 ! 180 ! Date of final time step 181 CALL calc_date( nitend, ditend_date ) 182 183 ! Background time for Jb referenced to ndate0 184 CALL calc_date( nitbkg_r, ditbkg_date ) 185 186 ! Background time for DI referenced to ndate0 187 CALL calc_date( nitdin_r, ditdin_date ) 188 189 ! IAU start time referenced to ndate0 190 CALL calc_date( nitiaustr_r, ditiaustr_date ) 191 192 ! IAU end time referenced to ndate0 193 CALL calc_date( nitiaufin_r, ditiaufin_date ) 194 186 195 IF(lwp) THEN 187 196 WRITE(numout,*) … … 198 207 WRITE(numout,*) ' ndastp = ', ndastp 199 208 WRITE(numout,*) ' ndate0 = ', ndate0 200 WRITE(numout,*) ' iitend_date = ', iitend_date 201 WRITE(numout,*) ' iitbkg_date = ', iitbkg_date 202 WRITE(numout,*) ' iitdin_date = ', iitdin_date 203 WRITE(numout,*) ' iitiaustr_date = ', iitiaustr_date 204 WRITE(numout,*) ' iitiaufin_date = ', iitiaufin_date 209 WRITE(numout,*) ' nn_time0 = ', nn_time0 210 WRITE(numout,*) ' ditend_date = ', ditend_date 211 WRITE(numout,*) ' ditbkg_date = ', ditbkg_date 212 WRITE(numout,*) ' ditdin_date = ', ditdin_date 213 WRITE(numout,*) ' ditiaustr_date = ', ditiaustr_date 214 WRITE(numout,*) ' ditiaufin_date = ', ditiaufin_date 205 215 ENDIF 206 216 207 IF ( nacc /= 0 ) &208 & CALL ctl_stop( ' nacc /= 0 and key_asminc :', &209 & ' Assimilation increments have only been implemented', &210 & ' for synchronous time stepping' )211 217 212 218 IF ( ( ln_asmdin ).AND.( ln_asmiau ) ) & … … 360 366 WRITE(numout,*) 361 367 WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', & 362 & ' between dates ', NINT( z_inc_dateb ),' and ', &363 & NINT( z_inc_datef )368 & ' between dates ', z_inc_dateb,' and ', & 369 & z_inc_datef 364 370 WRITE(numout,*) '~~~~~~~~~~~~' 365 371 ENDIF 366 372 367 IF ( ( NINT( z_inc_dateb ) < ndastp) &368 & .OR.( NINT( z_inc_datef ) > iitend_date ) ) &373 IF ( ( z_inc_dateb < ndastp + nn_time0*0.0001_wp ) & 374 & .OR.( z_inc_datef > ditend_date ) ) & 369 375 & CALL ctl_warn( ' Validity time of assimilation increments is ', & 370 376 & ' outside the assimilation interval' ) 371 377 372 IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &378 IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_date ) ) & 373 379 & CALL ctl_warn( ' Validity time of assimilation increments does ', & 374 380 & ' not agree with Direct Initialization time' ) … … 430 436 DO jt = 1, nn_divdmp 431 437 ! 432 DO jk = 1, jpkm1 438 DO jk = 1, jpkm1 ! hdiv = e1e1 * div 433 439 hdiv(:,:) = 0._wp 434 440 DO jj = 2, jpjm1 435 441 DO ji = fs_2, fs_jpim1 ! vector opt. 436 hdiv(ji,jj) = & 437 ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) & 438 - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) & 439 + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) & 440 - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) & 441 / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 442 hdiv(ji,jj) = ( e2u(ji ,jj) * e3u_n(ji ,jj,jk) * u_bkginc(ji ,jj,jk) & 443 & - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk) & 444 & + e1v(ji,jj ) * e3v_n(ji,jj ,jk) * v_bkginc(ji,jj ,jk) & 445 & - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk) ) / e3t_n(ji,jj,jk) 442 446 END DO 443 447 END DO … … 446 450 DO jj = 2, jpjm1 447 451 DO ji = fs_2, fs_jpim1 ! vector opt. 448 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1e2t(ji+1,jj) * hdiv(ji+1,jj) & 449 & - e1e2t(ji ,jj) * hdiv(ji ,jj) ) & 450 & * r1_e1u(ji,jj) * umask(ji,jj,jk) 451 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1e2t(ji,jj+1) * hdiv(ji,jj+1) & 452 & - e1e2t(ji,jj ) * hdiv(ji,jj ) ) & 453 & * r1_e2v(ji,jj) * vmask(ji,jj,jk) 452 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) & 453 & + 0.2_wp * ( hdiv(ji+1,jj) - hdiv(ji ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 454 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) & 455 & + 0.2_wp * ( hdiv(ji,jj+1) - hdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 454 456 END DO 455 457 END DO … … 490 492 IF(lwp) THEN 491 493 WRITE(numout,*) 492 WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', NINT( zdate_bkg ) 494 WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', & 495 & zdate_bkg 493 496 WRITE(numout,*) '~~~~~~~~~~~~' 494 497 ENDIF 495 498 ! 496 IF ( NINT( zdate_bkg ) /= iitdin_date ) &499 IF ( zdate_bkg /= ditdin_date ) & 497 500 & CALL ctl_warn( ' Validity time of assimilation background state does', & 498 501 & ' not agree with Direct Initialization time' ) … … 522 525 ! 523 526 END SUBROUTINE asm_inc_init 524 525 526 SUBROUTINE calc_date( kit000, kt, kdate0, kdate )527 !!----------------------------------------------------------------------528 !! *** ROUTINE calc_date ***529 !!530 !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.531 !!532 !! ** Method : Compute the calendar date YYYYMMDD at a given time step.533 !!534 !! ** Action :535 !!----------------------------------------------------------------------536 INTEGER, INTENT(IN) :: kit000 ! Initial time step537 INTEGER, INTENT(IN) :: kt ! Current time step referenced to kit000538 INTEGER, INTENT(IN) :: kdate0 ! Initial date539 INTEGER, INTENT(OUT) :: kdate ! Current date reference to kdate0540 !541 INTEGER :: iyea0 ! Initial year542 INTEGER :: imon0 ! Initial month543 INTEGER :: iday0 ! Initial day544 INTEGER :: iyea ! Current year545 INTEGER :: imon ! Current month546 INTEGER :: iday ! Current day547 INTEGER :: idaystp ! Number of days between initial and current date548 INTEGER :: idaycnt ! Day counter549 550 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year551 552 !-----------------------------------------------------------------------553 ! Compute the calendar date YYYYMMDD554 !-----------------------------------------------------------------------555 556 ! Initial date557 iyea0 = kdate0 / 10000558 imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100559 iday0 = kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )560 561 ! Check that kt >= kit000 - 1562 IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')563 564 ! If kt = kit000 - 1 then set the date to the restart date565 IF ( kt == kit000 - 1 ) THEN566 kdate = ndastp567 RETURN568 ENDIF569 570 ! Compute the number of days from the initial date571 idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )572 573 iday = iday0574 imon = imon0575 iyea = iyea0576 idaycnt = 0577 578 CALL calc_month_len( iyea, imonth_len )579 580 DO WHILE ( idaycnt < idaystp )581 iday = iday + 1582 IF ( iday > imonth_len(imon) ) THEN583 iday = 1584 imon = imon + 1585 ENDIF586 IF ( imon > 12 ) THEN587 imon = 1588 iyea = iyea + 1589 CALL calc_month_len( iyea, imonth_len ) ! update month lengths590 ENDIF591 idaycnt = idaycnt + 1592 END DO593 !594 kdate = iyea * 10000 + imon * 100 + iday595 !596 END SUBROUTINE597 598 599 SUBROUTINE calc_month_len( iyear, imonth_len )600 !!----------------------------------------------------------------------601 !! *** ROUTINE calc_month_len ***602 !!603 !! ** Purpose : Compute the number of days in a months given a year.604 !!605 !! ** Method :606 !!----------------------------------------------------------------------607 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year608 INTEGER :: iyear !: year609 !!----------------------------------------------------------------------610 !611 ! length of the month of the current year (from nleapy, read in namelist)612 IF ( nleapy < 2 ) THEN613 imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)614 IF ( nleapy == 1 ) THEN ! we are using calendar with leap years615 IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN616 imonth_len(2) = 29617 ENDIF618 ENDIF619 ELSE620 imonth_len(:) = nleapy ! all months with nleapy days per year621 ENDIF622 !623 END SUBROUTINE624 625 626 527 SUBROUTINE tra_asm_inc( kt ) 627 528 !!---------------------------------------------------------------------- … … 645 546 ! used to prevent the applied increments taking the temperature below the local freezing point 646 547 DO jk = 1, jpkm1 647 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) )548 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) ) 648 549 END DO 649 550 ! … … 726 627 !!gm 727 628 728 IF( ln_zps .AND. .NOT. lk_c1d ) THEN ! Partial steps: before horizontal gradient 729 IF(ln_isfcav) THEN ! ocean cavities: top and bottom cells (ISF) 730 CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi, & 731 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 732 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 733 ELSE ! no ocean cavities: bottom cells 734 CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! 735 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 736 ENDIF 737 ENDIF 738 ! 629 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 630 & CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 631 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 632 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 633 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 634 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level 635 739 636 DEALLOCATE( t_bkginc ) 740 637 DEALLOCATE( s_bkginc ) … … 874 771 ! 875 772 sshb(:,:) = sshn(:,:) ! Update before fields 876 ! 877 IF( lk_vvl ) THEN 878 DO jk = 1, jpk 879 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 880 END DO 881 ENDIF 773 e3t_b(:,:,:) = e3t_n(:,:,:) 774 !!gm why not e3u_b, e3v_b, gdept_b ???? 882 775 ! 883 776 DEALLOCATE( ssh_bkg )
Note: See TracChangeset
for help on using the changeset viewer.
