- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r5541 r6808 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 !!----------------------------------------------------------------------16 !! 'key_asminc' : Switch on the assimilation increment interface17 17 !!---------------------------------------------------------------------- 18 18 !! asm_inc_init : Initialize the increment arrays and IAU weights 19 !! calc_date : Compute the calendar date YYYYMMDD on a given step20 19 !! tra_asm_inc : Apply the tracer (T and S) increments 21 20 !! dyn_asm_inc : Apply the dynamic (u and v) increments … … 28 27 USE domvvl ! domain: variable volume level 29 28 USE oce ! Dynamics and active tracers defined in memory 30 USE ldfdyn _oce ! ocean dynamics: lateral physics29 USE ldfdyn ! lateral diffusion: eddy viscosity coefficients 31 30 USE eosbn2 ! Equation of state - in situ and potential density 32 31 USE zpshde ! Partial step : Horizontal Derivative … … 40 39 #endif 41 40 USE sbc_oce ! Surface boundary condition variables. 41 USE diaobs, ONLY: calc_date ! Compute the calendar date on a given step 42 42 43 43 IMPLICIT NONE … … 45 45 46 46 PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights 47 PUBLIC calc_date !: Compute the calendar date YYYYMMDD on a given step48 47 PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments 49 48 PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments … … 56 55 LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE. !: No assimilation increments 57 56 #endif 58 LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.!: No output of the background state fields59 LOGICAL, PUBLIC :: ln_asmiau = .FALSE.!: No applying forcing with an assimilation increment60 LOGICAL, PUBLIC :: ln_asmdin = .FALSE.!: No direct initialization61 LOGICAL, PUBLIC :: ln_trainc = .FALSE.!: No tracer (T and S) assimilation increments62 LOGICAL, PUBLIC :: ln_dyninc = .FALSE.!: No dynamics (u and v) assimilation increments63 LOGICAL, PUBLIC :: ln_sshinc = .FALSE.!: No sea surface height assimilation increment64 LOGICAL, PUBLIC :: ln_seaiceinc 65 LOGICAL, PUBLIC :: ln_salfix = .FALSE.!: Apply minimum salinity check57 LOGICAL, PUBLIC :: ln_bkgwri !: No output of the background state fields 58 LOGICAL, PUBLIC :: ln_asmiau !: No applying forcing with an assimilation increment 59 LOGICAL, PUBLIC :: ln_asmdin !: No direct initialization 60 LOGICAL, PUBLIC :: ln_trainc !: No tracer (T and S) assimilation increments 61 LOGICAL, PUBLIC :: ln_dyninc !: No dynamics (u and v) assimilation increments 62 LOGICAL, PUBLIC :: ln_sshinc !: No sea surface height assimilation increment 63 LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment 64 LOGICAL, PUBLIC :: ln_salfix !: Apply minimum salinity check 66 65 LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 67 INTEGER, PUBLIC :: nn_divdmp 66 INTEGER, PUBLIC :: nn_divdmp !: Apply divergence damping filter nn_divdmp times 68 67 69 68 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkg , s_bkg !: Background temperature and salinity … … 89 88 90 89 !! * Substitutions 91 # include "domzgr_substitute.h90"92 # include "ldfdyn_substitute.h90"93 90 # include "vectopt_loop_substitute.h90" 94 91 !!---------------------------------------------------------------------- 95 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)92 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 96 93 !! $Id$ 97 94 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 114 111 INTEGER :: iiauper ! Number of time steps in the IAU period 115 112 INTEGER :: icycper ! Number of time steps in the cycle 116 INTEGER :: iitend_date ! Date YYYYMMDDof final time step117 INTEGER :: iitbkg_date ! Date YYYYMMDDof background time step for Jb term118 INTEGER :: iitdin_date ! Date YYYYMMDDof background time step for DI119 INTEGER :: iitiaustr_date ! Date YYYYMMDDof IAU interval start time step120 INTEGER :: iitiaufin_date ! Date YYYYMMDDof IAU interval final time step121 ! 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 122 119 REAL(wp) :: znorm ! Normalization factor for IAU weights 123 120 REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights (should be equal to one) … … 139 136 ! Read Namelist nam_asminc : assimilation increment interface 140 137 !----------------------------------------------------------------------- 141 ln_seaiceinc = .FALSE.138 ln_seaiceinc = .FALSE. 142 139 ln_temnofreeze = .FALSE. 143 140 … … 181 178 icycper = nitend - nit000 + 1 ! Cycle interval length 182 179 183 CALL calc_date( nit000, nitend , ndate0, iitend_date ) ! Date of final time step 184 CALL calc_date( nit000, nitbkg_r , ndate0, iitbkg_date ) ! Background time for Jb referenced to ndate0 185 CALL calc_date( nit000, nitdin_r , ndate0, iitdin_date ) ! Background time for DI referenced to ndate0 186 CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date ) ! IAU start time referenced to ndate0 187 CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date ) ! IAU end time referenced to ndate0 188 ! 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 189 195 IF(lwp) THEN 190 196 WRITE(numout,*) … … 201 207 WRITE(numout,*) ' ndastp = ', ndastp 202 208 WRITE(numout,*) ' ndate0 = ', ndate0 203 WRITE(numout,*) ' iitend_date = ', iitend_date 204 WRITE(numout,*) ' iitbkg_date = ', iitbkg_date 205 WRITE(numout,*) ' iitdin_date = ', iitdin_date 206 WRITE(numout,*) ' iitiaustr_date = ', iitiaustr_date 207 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 208 215 ENDIF 209 216 210 IF ( nacc /= 0 ) &211 & CALL ctl_stop( ' nacc /= 0 and key_asminc :', &212 & ' Assimilation increments have only been implemented', &213 & ' for synchronous time stepping' )214 217 215 218 IF ( ( ln_asmdin ).AND.( ln_asmiau ) ) & … … 363 366 WRITE(numout,*) 364 367 WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', & 365 & ' between dates ', NINT( z_inc_dateb ),' and ', &366 & NINT( z_inc_datef )368 & ' between dates ', z_inc_dateb,' and ', & 369 & z_inc_datef 367 370 WRITE(numout,*) '~~~~~~~~~~~~' 368 371 ENDIF 369 372 370 IF ( ( NINT( z_inc_dateb ) < ndastp) &371 & .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 ) ) & 372 375 & CALL ctl_warn( ' Validity time of assimilation increments is ', & 373 376 & ' outside the assimilation interval' ) 374 377 375 IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &378 IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_date ) ) & 376 379 & CALL ctl_warn( ' Validity time of assimilation increments does ', & 377 380 & ' not agree with Direct Initialization time' ) … … 428 431 429 432 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 430 431 CALL wrk_alloc(jpi,jpj,hdiv) 432 433 DO jt = 1, nn_divdmp 434 435 DO jk = 1, jpkm1 436 433 ! 434 CALL wrk_alloc( jpi,jpj, hdiv ) 435 ! 436 DO jt = 1, nn_divdmp 437 ! 438 DO jk = 1, jpkm1 ! hdiv = e1e1 * div 437 439 hdiv(:,:) = 0._wp 438 439 440 DO jj = 2, jpjm1 440 441 DO ji = fs_2, fs_jpim1 ! vector opt. 441 hdiv(ji,jj) = & 442 ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) & 443 - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) & 444 + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) & 445 - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) & 446 / ( e1t(ji,jj) * e2t(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) 447 446 END DO 448 447 END DO 449 450 448 CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) 451 449 ! 452 450 DO jj = 2, jpjm1 453 451 DO ji = fs_2, fs_jpim1 ! vector opt. 454 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj) & 455 - e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) & 456 / e1u(ji,jj) * umask(ji,jj,jk) 457 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1) & 458 - e1t(ji,jj )*e2t(ji,jj ) * hdiv(ji,jj ) ) & 459 / 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) 460 456 END DO 461 457 END DO 462 463 458 END DO 464 459 ! 465 460 END DO 466 467 CALL wrk_dealloc( jpi,jpj,hdiv)468 461 ! 462 CALL wrk_dealloc( jpi,jpj, hdiv ) 463 ! 469 464 ENDIF 470 471 472 465 473 466 !----------------------------------------------------------------------- … … 476 469 477 470 IF ( ln_asmdin ) THEN 478 471 ! 479 472 ALLOCATE( t_bkg(jpi,jpj,jpk) ) 480 473 ALLOCATE( s_bkg(jpi,jpj,jpk) ) … … 482 475 ALLOCATE( v_bkg(jpi,jpj,jpk) ) 483 476 ALLOCATE( ssh_bkg(jpi,jpj) ) 484 485 t_bkg(:,:,:) = 0. 0486 s_bkg(:,:,:) = 0. 0487 u_bkg(:,:,:) = 0. 0488 v_bkg(:,:,:) = 0. 0489 ssh_bkg(:,:) = 0. 0490 477 ! 478 t_bkg(:,:,:) = 0._wp 479 s_bkg(:,:,:) = 0._wp 480 u_bkg(:,:,:) = 0._wp 481 v_bkg(:,:,:) = 0._wp 482 ssh_bkg(:,:) = 0._wp 483 ! 491 484 !-------------------------------------------------------------------- 492 485 ! Read from file the background state at analysis time 493 486 !-------------------------------------------------------------------- 494 487 ! 495 488 CALL iom_open( c_asmdin, inum ) 496 489 ! 497 490 CALL iom_get( inum, 'rdastp', zdate_bkg ) 498 491 ! 499 492 IF(lwp) THEN 500 493 WRITE(numout,*) 501 494 WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', & 502 & NINT( zdate_bkg )495 & zdate_bkg 503 496 WRITE(numout,*) '~~~~~~~~~~~~' 504 497 ENDIF 505 506 IF ( NINT( zdate_bkg ) /= iitdin_date ) &498 ! 499 IF ( zdate_bkg /= ditdin_date ) & 507 500 & CALL ctl_warn( ' Validity time of assimilation background state does', & 508 501 & ' not agree with Direct Initialization time' ) 509 502 ! 510 503 IF ( ln_trainc ) THEN 511 504 CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg ) … … 514 507 s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:) 515 508 ENDIF 516 509 ! 517 510 IF ( ln_dyninc ) THEN 518 511 CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg ) … … 521 514 v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:) 522 515 ENDIF 523 516 ! 524 517 IF ( ln_sshinc ) THEN 525 518 CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg ) 526 519 ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1) 527 520 ENDIF 528 521 ! 529 522 CALL iom_close( inum ) 530 523 ! 531 524 ENDIF 532 525 ! 533 526 END SUBROUTINE asm_inc_init 534 535 536 SUBROUTINE calc_date( kit000, kt, kdate0, kdate )537 !!----------------------------------------------------------------------538 !! *** ROUTINE calc_date ***539 !!540 !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.541 !!542 !! ** Method : Compute the calendar date YYYYMMDD at a given time step.543 !!544 !! ** Action :545 !!----------------------------------------------------------------------546 INTEGER, INTENT(IN) :: kit000 ! Initial time step547 INTEGER, INTENT(IN) :: kt ! Current time step referenced to kit000548 INTEGER, INTENT(IN) :: kdate0 ! Initial date549 INTEGER, INTENT(OUT) :: kdate ! Current date reference to kdate0550 !551 INTEGER :: iyea0 ! Initial year552 INTEGER :: imon0 ! Initial month553 INTEGER :: iday0 ! Initial day554 INTEGER :: iyea ! Current year555 INTEGER :: imon ! Current month556 INTEGER :: iday ! Current day557 INTEGER :: idaystp ! Number of days between initial and current date558 INTEGER :: idaycnt ! Day counter559 560 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year561 562 !-----------------------------------------------------------------------563 ! Compute the calendar date YYYYMMDD564 !-----------------------------------------------------------------------565 566 ! Initial date567 iyea0 = kdate0 / 10000568 imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100569 iday0 = kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )570 571 ! Check that kt >= kit000 - 1572 IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')573 574 ! If kt = kit000 - 1 then set the date to the restart date575 IF ( kt == kit000 - 1 ) THEN576 577 kdate = ndastp578 RETURN579 580 ENDIF581 582 ! Compute the number of days from the initial date583 idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )584 585 iday = iday0586 imon = imon0587 iyea = iyea0588 idaycnt = 0589 590 CALL calc_month_len( iyea, imonth_len )591 592 DO WHILE ( idaycnt < idaystp )593 iday = iday + 1594 IF ( iday > imonth_len(imon) ) THEN595 iday = 1596 imon = imon + 1597 ENDIF598 IF ( imon > 12 ) THEN599 imon = 1600 iyea = iyea + 1601 CALL calc_month_len( iyea, imonth_len ) ! update month lengths602 ENDIF603 idaycnt = idaycnt + 1604 END DO605 !606 kdate = iyea * 10000 + imon * 100 + iday607 !608 END SUBROUTINE609 610 611 SUBROUTINE calc_month_len( iyear, imonth_len )612 !!----------------------------------------------------------------------613 !! *** ROUTINE calc_month_len ***614 !!615 !! ** Purpose : Compute the number of days in a months given a year.616 !!617 !! ** Method :618 !!----------------------------------------------------------------------619 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year620 INTEGER :: iyear !: year621 !!----------------------------------------------------------------------622 !623 ! length of the month of the current year (from nleapy, read in namelist)624 IF ( nleapy < 2 ) THEN625 imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)626 IF ( nleapy == 1 ) THEN ! we are using calendar with leap years627 IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN628 imonth_len(2) = 29629 ENDIF630 ENDIF631 ELSE632 imonth_len(:) = nleapy ! all months with nleapy days per year633 ENDIF634 !635 END SUBROUTINE636 637 638 527 SUBROUTINE tra_asm_inc( kt ) 639 528 !!---------------------------------------------------------------------- … … 646 535 !! ** Action : 647 536 !!---------------------------------------------------------------------- 648 INTEGER, INTENT(IN) :: kt! Current time step649 ! 650 INTEGER :: ji,jj,jk651 INTEGER :: it537 INTEGER, INTENT(IN) :: kt ! Current time step 538 ! 539 INTEGER :: ji, jj, jk 540 INTEGER :: it 652 541 REAL(wp) :: zincwgt ! IAU weight for current time step 653 542 REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 654 543 !!---------------------------------------------------------------------- 655 544 ! 656 545 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 657 546 ! used to prevent the applied increments taking the temperature below the local freezing point 658 659 547 DO jk = 1, jpkm1 660 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) )548 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) ) 661 549 END DO 662 663 IF ( ln_asmiau ) THEN 664 665 !-------------------------------------------------------------------- 666 ! Incremental Analysis Updating 667 !-------------------------------------------------------------------- 668 550 ! 551 ! !-------------------------------------- 552 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 553 ! !-------------------------------------- 554 ! 669 555 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 670 556 ! 671 557 it = kt - nit000 + 1 672 558 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 673 559 ! 674 560 IF(lwp) THEN 675 561 WRITE(numout,*) … … 677 563 WRITE(numout,*) '~~~~~~~~~~~~' 678 564 ENDIF 679 565 ! 680 566 ! Update the tracer tendencies 681 567 DO jk = 1, jpkm1 … … 700 586 ENDIF 701 587 END DO 702 703 ENDIF 704 588 ! 589 ENDIF 590 ! 705 591 IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work 706 592 DEALLOCATE( t_bkginc ) 707 593 DEALLOCATE( s_bkginc ) 708 594 ENDIF 709 710 711 ELSEIF ( ln_asmdin ) THEN 712 713 !-------------------------------------------------------------------- 714 ! Direct Initialization 715 !-------------------------------------------------------------------- 716 595 ! !-------------------------------------- 596 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 597 ! !-------------------------------------- 598 ! 717 599 IF ( kt == nitdin_r ) THEN 718 600 ! 719 601 neuler = 0 ! Force Euler forward step 720 602 ! 721 603 ! Initialize the now fields with the background + increment 722 604 IF (ln_temnofreeze) THEN … … 745 627 !!gm 746 628 747 748 629 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 749 630 & CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 750 631 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 751 632 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 752 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF) 753 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 754 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 755 756 #if defined key_zdfkpp 757 CALL eos( tsn, rhd, fsdept_n(:,:,:) ) ! Compute rhd 758 !!gm fabien CALL eos( tsn, rhd ) ! Compute rhd 759 #endif 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 760 635 761 636 DEALLOCATE( t_bkginc ) … … 767 642 ENDIF 768 643 ! Perhaps the following call should be in step 769 IF 644 IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment 770 645 ! 771 646 END SUBROUTINE tra_asm_inc … … 788 663 REAL(wp) :: zincwgt ! IAU weight for current time step 789 664 !!---------------------------------------------------------------------- 790 791 IF ( ln_asmiau ) THEN 792 793 !-------------------------------------------------------------------- 794 ! Incremental Analysis Updating 795 !-------------------------------------------------------------------- 796 665 ! 666 ! !-------------------------------------------- 667 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 668 ! !-------------------------------------------- 669 ! 797 670 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 798 671 ! 799 672 it = kt - nit000 + 1 800 673 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 801 674 ! 802 675 IF(lwp) THEN 803 676 WRITE(numout,*) 804 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', & 805 & kt,' with IAU weight = ', wgtiau(it) 677 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 806 678 WRITE(numout,*) '~~~~~~~~~~~~' 807 679 ENDIF 808 680 ! 809 681 ! Update the dynamic tendencies 810 682 DO jk = 1, jpkm1 … … 812 684 va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt 813 685 END DO 814 686 ! 815 687 IF ( kt == nitiaufin_r ) THEN 816 688 DEALLOCATE( u_bkginc ) 817 689 DEALLOCATE( v_bkginc ) 818 690 ENDIF 819 820 ENDIF 821 822 ELSEIF ( ln_asmdin ) THEN 823 824 !-------------------------------------------------------------------- 825 ! Direct Initialization 826 !-------------------------------------------------------------------- 827 691 ! 692 ENDIF 693 ! !----------------------------------------- 694 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 695 ! !----------------------------------------- 696 ! 828 697 IF ( kt == nitdin_r ) THEN 829 698 ! 830 699 neuler = 0 ! Force Euler forward step 831 700 ! 832 701 ! Initialize the now fields with the background + increment 833 702 un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:) 834 703 vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 835 704 ! 836 705 ub(:,:,:) = un(:,:,:) ! Update before fields 837 706 vb(:,:,:) = vn(:,:,:) 838 707 ! 839 708 DEALLOCATE( u_bkg ) 840 709 DEALLOCATE( v_bkg ) … … 864 733 REAL(wp) :: zincwgt ! IAU weight for current time step 865 734 !!---------------------------------------------------------------------- 866 867 IF ( ln_asmiau ) THEN 868 869 !-------------------------------------------------------------------- 870 ! Incremental Analysis Updating 871 !-------------------------------------------------------------------- 872 735 ! 736 ! !----------------------------------------- 737 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 738 ! !----------------------------------------- 739 ! 873 740 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 874 741 ! 875 742 it = kt - nit000 + 1 876 743 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 877 744 ! 878 745 IF(lwp) THEN 879 746 WRITE(numout,*) … … 882 749 WRITE(numout,*) '~~~~~~~~~~~~' 883 750 ENDIF 884 751 ! 885 752 ! Save the tendency associated with the IAU weighted SSH increment 886 753 ! (applied in dynspg.*) … … 891 758 DEALLOCATE( ssh_bkginc ) 892 759 ENDIF 893 894 ENDIF 895 896 ELSEIF ( ln_asmdin ) THEN 897 898 !-------------------------------------------------------------------- 899 ! Direct Initialization 900 !-------------------------------------------------------------------- 901 760 ! 761 ENDIF 762 ! !----------------------------------------- 763 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 764 ! !----------------------------------------- 765 ! 902 766 IF ( kt == nitdin_r ) THEN 903 904 neuler = 0 ! Force Euler forward step 905 906 ! Initialize the now fields the background + increment 907 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 908 909 ! Update before fields 910 sshb(:,:) = sshn(:,:) 911 912 IF( lk_vvl ) THEN 913 DO jk = 1, jpk 914 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 915 END DO 916 ENDIF 917 767 ! 768 neuler = 0 ! Force Euler forward step 769 ! 770 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment 771 ! 772 sshb(:,:) = sshn(:,:) ! Update before fields 773 e3t_b(:,:,:) = e3t_n(:,:,:) 774 !!gm why not e3u_b, e3v_b, gdept_b ???? 775 ! 918 776 DEALLOCATE( ssh_bkg ) 919 777 DEALLOCATE( ssh_bkginc ) 920 778 ! 921 779 ENDIF 922 780 ! … … 937 795 !! 938 796 !!---------------------------------------------------------------------- 939 IMPLICIT NONE 940 ! 941 INTEGER, INTENT(in) :: kt ! Current time step 797 INTEGER, INTENT(in) :: kt ! Current time step 942 798 INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation 943 799 ! … … 949 805 #endif 950 806 !!---------------------------------------------------------------------- 951 952 IF ( ln_asmiau ) THEN 953 954 !-------------------------------------------------------------------- 955 ! Incremental Analysis Updating 956 !-------------------------------------------------------------------- 957 807 ! 808 ! !----------------------------------------- 809 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 810 ! !----------------------------------------- 811 ! 958 812 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 959 813 ! 960 814 it = kt - nit000 + 1 961 815 zincwgt = wgtiau(it) ! IAU weight for the current time step 962 816 ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 963 817 ! 964 818 IF(lwp) THEN 965 819 WRITE(numout,*) 966 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 967 & kt,' with IAU weight = ', wgtiau(it) 820 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 968 821 WRITE(numout,*) '~~~~~~~~~~~~' 969 822 ENDIF 970 823 ! 971 824 ! Sea-ice : LIM-3 case (to add) 972 825 ! 973 826 #if defined key_lim2 974 827 ! Sea-ice : LIM-2 case … … 1008 861 1009 862 #if defined key_cice && defined key_asminc 1010 ! Sea-ice : CICE case. Zero ice increment tendency into CICE 1011 ndaice_da(:,:) = 0.0_wp 1012 #endif 1013 1014 ENDIF 1015 1016 ELSEIF ( ln_asmdin ) THEN 1017 1018 !-------------------------------------------------------------------- 1019 ! Direct Initialization 1020 !-------------------------------------------------------------------- 1021 863 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 864 #endif 865 866 ENDIF 867 ! !----------------------------------------- 868 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 869 ! !----------------------------------------- 870 ! 1022 871 IF ( kt == nitdin_r ) THEN 1023 872 ! 1024 873 neuler = 0 ! Force Euler forward step 1025 874 ! 1026 875 ! Sea-ice : LIM-3 case (to add) 1027 876 ! 1028 877 #if defined key_lim2 1029 878 ! Sea-ice : LIM-2 case. … … 1041 890 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 1042 891 ELSEWHERE 1043 zhicifinc(:,:) = 0. 0_wp892 zhicifinc(:,:) = 0._wp 1044 893 END WHERE 1045 894 ! … … 1050 899 ! seaice salinity balancing (to add) 1051 900 #endif 1052 901 ! 1053 902 #if defined key_cice && defined key_asminc 1054 903 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 1055 904 ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 1056 905 #endif 1057 IF ( .NOT. PRESENT(kindic) ) THEN1058 DEALLOCATE( seaice_bkginc )1059 END IF1060 906 IF ( .NOT. PRESENT(kindic) ) THEN 907 DEALLOCATE( seaice_bkginc ) 908 END IF 909 ! 1061 910 ELSE 1062 911 ! 1063 912 #if defined key_cice && defined key_asminc 1064 ! Sea-ice : CICE case. Zero ice increment tendency into CICE1065 ndaice_da(:,:) = 0.0_wp 1066 #endif 1067 913 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 914 915 #endif 916 ! 1068 917 ENDIF 1069 918 … … 1142 991 ! 1143 992 !#endif 1144 993 ! 1145 994 ENDIF 1146 995 ! 1147 996 END SUBROUTINE seaice_asm_inc 1148 997
Note: See TracChangeset
for help on using the changeset viewer.
