- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r4624 r6225 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 … … 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) … … 109 106 !! ** Action : 110 107 !!---------------------------------------------------------------------- 111 INTEGER :: ji, jj, jk 112 INTEGER :: jt 113 INTEGER :: imid 114 INTEGER :: inum 108 INTEGER :: ji, jj, jk, jt ! dummy loop indices 109 INTEGER :: imid, inum ! local integers 110 INTEGER :: ios ! Local integer output status for namelist read 115 111 INTEGER :: iiauper ! Number of time steps in the IAU period 116 112 INTEGER :: icycper ! Number of time steps in the cycle 117 INTEGER :: iitend_date ! Date YYYYMMDD of final time step 118 INTEGER :: iitbkg_date ! Date YYYYMMDD of background time step for Jb term 119 INTEGER :: iitdin_date ! Date YYYYMMDD of background time step for DI 120 INTEGER :: iitiaustr_date ! Date YYYYMMDD of IAU interval start time step 121 INTEGER :: iitiaufin_date ! Date YYYYMMDD of IAU interval final time step 122 INTEGER :: ios ! Local integer output status for namelist read 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 123 118 124 119 REAL(wp) :: znorm ! Normalization factor for IAU weights 125 REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights 126 ! (should be equal to one) 120 REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights (should be equal to one) 127 121 REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid 128 122 REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid 129 123 REAL(wp) :: zdate_bkg ! Date in background state file for DI 130 124 REAL(wp) :: zdate_inc ! Time axis in increments file 131 132 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv125 ! 126 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv ! 2D workspace 133 127 !! 134 128 NAMELIST/nam_asminc/ ln_bkgwri, & … … 136 130 & ln_asmdin, ln_asmiau, & 137 131 & nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, & 138 & ln_salfix, salfixmin, & 139 & nn_divdmp 132 & ln_salfix, salfixmin, nn_divdmp 140 133 !!---------------------------------------------------------------------- 141 134 … … 143 136 ! Read Namelist nam_asminc : assimilation increment interface 144 137 !----------------------------------------------------------------------- 145 146 ln_seaiceinc = .FALSE. 138 ln_seaiceinc = .FALSE. 147 139 ln_temnofreeze = .FALSE. 148 140 … … 187 179 188 180 ! Date of final time step 189 CALL calc_date( nit 000, nitend, ndate0, iitend_date )181 CALL calc_date( nitend, ditend_date ) 190 182 191 183 ! Background time for Jb referenced to ndate0 192 CALL calc_date( nit 000, nitbkg_r, ndate0, iitbkg_date )184 CALL calc_date( nitbkg_r, ditbkg_date ) 193 185 194 186 ! Background time for DI referenced to ndate0 195 CALL calc_date( nit 000, nitdin_r, ndate0, iitdin_date )187 CALL calc_date( nitdin_r, ditdin_date ) 196 188 197 189 ! IAU start time referenced to ndate0 198 CALL calc_date( nit 000, nitiaustr_r, ndate0, iitiaustr_date )190 CALL calc_date( nitiaustr_r, ditiaustr_date ) 199 191 200 192 ! IAU end time referenced to ndate0 201 CALL calc_date( nit 000, nitiaufin_r, ndate0, iitiaufin_date )193 CALL calc_date( nitiaufin_r, ditiaufin_date ) 202 194 203 195 IF(lwp) THEN … … 215 207 WRITE(numout,*) ' ndastp = ', ndastp 216 208 WRITE(numout,*) ' ndate0 = ', ndate0 217 WRITE(numout,*) ' iitend_date = ', iitend_date 218 WRITE(numout,*) ' iitbkg_date = ', iitbkg_date 219 WRITE(numout,*) ' iitdin_date = ', iitdin_date 220 WRITE(numout,*) ' iitiaustr_date = ', iitiaustr_date 221 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 222 215 ENDIF 223 216 224 IF ( nacc /= 0 ) &225 & CALL ctl_stop( ' nacc /= 0 and key_asminc :', &226 & ' Assimilation increments have only been implemented', &227 & ' for synchronous time stepping' )228 217 229 218 IF ( ( ln_asmdin ).AND.( ln_asmiau ) ) & … … 236 225 & ' but ln_asmdin and ln_asmiau are both set to .false. :', & 237 226 & ' Inconsistent options') 238 239 IF ( ( ln_bkgwri ).AND.( ( ln_asmdin ).OR.( ln_asmiau ) ) ) &240 & CALL ctl_stop( ' ln_bkgwri and either ln_asmdin or ln_asmiau are set to .true.:', &241 & ' The background state must be written before applying the increments')242 227 243 228 IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) & … … 381 366 WRITE(numout,*) 382 367 WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', & 383 & ' between dates ', NINT( z_inc_dateb ),' and ', &384 & NINT( z_inc_datef )368 & ' between dates ', z_inc_dateb,' and ', & 369 & z_inc_datef 385 370 WRITE(numout,*) '~~~~~~~~~~~~' 386 371 ENDIF 387 372 388 IF ( ( NINT( z_inc_dateb ) < ndastp) &389 & .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 ) ) & 390 375 & CALL ctl_warn( ' Validity time of assimilation increments is ', & 391 376 & ' outside the assimilation interval' ) 392 377 393 IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &378 IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_date ) ) & 394 379 & CALL ctl_warn( ' Validity time of assimilation increments does ', & 395 380 & ' not agree with Direct Initialization time' ) … … 446 431 447 432 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 448 449 CALL wrk_alloc(jpi,jpj,hdiv) 450 451 DO jt = 1, nn_divdmp 452 453 DO jk = 1, jpkm1 454 433 ! 434 CALL wrk_alloc( jpi,jpj, hdiv ) 435 ! 436 DO jt = 1, nn_divdmp 437 ! 438 DO jk = 1, jpkm1 ! hdiv = e1e1 * div 455 439 hdiv(:,:) = 0._wp 456 457 440 DO jj = 2, jpjm1 458 441 DO ji = fs_2, fs_jpim1 ! vector opt. 459 hdiv(ji,jj) = & 460 ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) & 461 - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) & 462 + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) & 463 - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) & 464 / ( 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) 465 446 END DO 466 447 END DO 467 468 448 CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) 469 449 ! 470 450 DO jj = 2, jpjm1 471 451 DO ji = fs_2, fs_jpim1 ! vector opt. 472 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) & 473 - e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) & 474 / e1u(ji,jj) * umask(ji,jj,jk) 475 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) & 476 - e1t(ji,jj )*e2t(ji,jj ) * hdiv(ji,jj ) ) & 477 / 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) 478 456 END DO 479 457 END DO 480 481 458 END DO 482 459 ! 483 460 END DO 484 485 CALL wrk_dealloc( jpi,jpj,hdiv)486 461 ! 462 CALL wrk_dealloc( jpi,jpj, hdiv ) 463 ! 487 464 ENDIF 488 489 490 465 491 466 !----------------------------------------------------------------------- … … 494 469 495 470 IF ( ln_asmdin ) THEN 496 471 ! 497 472 ALLOCATE( t_bkg(jpi,jpj,jpk) ) 498 473 ALLOCATE( s_bkg(jpi,jpj,jpk) ) … … 500 475 ALLOCATE( v_bkg(jpi,jpj,jpk) ) 501 476 ALLOCATE( ssh_bkg(jpi,jpj) ) 502 503 t_bkg(:,:,:) = 0. 0504 s_bkg(:,:,:) = 0. 0505 u_bkg(:,:,:) = 0. 0506 v_bkg(:,:,:) = 0. 0507 ssh_bkg(:,:) = 0. 0508 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 ! 509 484 !-------------------------------------------------------------------- 510 485 ! Read from file the background state at analysis time 511 486 !-------------------------------------------------------------------- 512 487 ! 513 488 CALL iom_open( c_asmdin, inum ) 514 489 ! 515 490 CALL iom_get( inum, 'rdastp', zdate_bkg ) 516 491 ! 517 492 IF(lwp) THEN 518 493 WRITE(numout,*) 519 494 WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', & 520 & NINT( zdate_bkg )495 & zdate_bkg 521 496 WRITE(numout,*) '~~~~~~~~~~~~' 522 497 ENDIF 523 524 IF ( NINT( zdate_bkg ) /= iitdin_date ) &498 ! 499 IF ( zdate_bkg /= ditdin_date ) & 525 500 & CALL ctl_warn( ' Validity time of assimilation background state does', & 526 501 & ' not agree with Direct Initialization time' ) 527 502 ! 528 503 IF ( ln_trainc ) THEN 529 504 CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg ) … … 532 507 s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:) 533 508 ENDIF 534 509 ! 535 510 IF ( ln_dyninc ) THEN 536 511 CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg ) … … 539 514 v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:) 540 515 ENDIF 541 516 ! 542 517 IF ( ln_sshinc ) THEN 543 518 CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg ) 544 519 ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1) 545 520 ENDIF 546 521 ! 547 522 CALL iom_close( inum ) 548 523 ! 549 524 ENDIF 550 525 ! 551 526 END SUBROUTINE asm_inc_init 552 553 554 SUBROUTINE calc_date( kit000, kt, kdate0, kdate )555 !!----------------------------------------------------------------------556 !! *** ROUTINE calc_date ***557 !!558 !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.559 !!560 !! ** Method : Compute the calendar date YYYYMMDD at a given time step.561 !!562 !! ** Action :563 !!----------------------------------------------------------------------564 INTEGER, INTENT(IN) :: kit000 ! Initial time step565 INTEGER, INTENT(IN) :: kt ! Current time step referenced to kit000566 INTEGER, INTENT(IN) :: kdate0 ! Initial date567 INTEGER, INTENT(OUT) :: kdate ! Current date reference to kdate0568 !569 INTEGER :: iyea0 ! Initial year570 INTEGER :: imon0 ! Initial month571 INTEGER :: iday0 ! Initial day572 INTEGER :: iyea ! Current year573 INTEGER :: imon ! Current month574 INTEGER :: iday ! Current day575 INTEGER :: idaystp ! Number of days between initial and current date576 INTEGER :: idaycnt ! Day counter577 578 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year579 580 !-----------------------------------------------------------------------581 ! Compute the calendar date YYYYMMDD582 !-----------------------------------------------------------------------583 584 ! Initial date585 iyea0 = kdate0 / 10000586 imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100587 iday0 = kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )588 589 ! Check that kt >= kit000 - 1590 IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')591 592 ! If kt = kit000 - 1 then set the date to the restart date593 IF ( kt == kit000 - 1 ) THEN594 595 kdate = ndastp596 RETURN597 598 ENDIF599 600 ! Compute the number of days from the initial date601 idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )602 603 iday = iday0604 imon = imon0605 iyea = iyea0606 idaycnt = 0607 608 CALL calc_month_len( iyea, imonth_len )609 610 DO WHILE ( idaycnt < idaystp )611 iday = iday + 1612 IF ( iday > imonth_len(imon) ) THEN613 iday = 1614 imon = imon + 1615 ENDIF616 IF ( imon > 12 ) THEN617 imon = 1618 iyea = iyea + 1619 CALL calc_month_len( iyea, imonth_len ) ! update month lengths620 ENDIF621 idaycnt = idaycnt + 1622 END DO623 !624 kdate = iyea * 10000 + imon * 100 + iday625 !626 END SUBROUTINE627 628 629 SUBROUTINE calc_month_len( iyear, imonth_len )630 !!----------------------------------------------------------------------631 !! *** ROUTINE calc_month_len ***632 !!633 !! ** Purpose : Compute the number of days in a months given a year.634 !!635 !! ** Method :636 !!----------------------------------------------------------------------637 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year638 INTEGER :: iyear !: year639 !!----------------------------------------------------------------------640 !641 ! length of the month of the current year (from nleapy, read in namelist)642 IF ( nleapy < 2 ) THEN643 imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)644 IF ( nleapy == 1 ) THEN ! we are using calendar with leap years645 IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN646 imonth_len(2) = 29647 ENDIF648 ENDIF649 ELSE650 imonth_len(:) = nleapy ! all months with nleapy days per year651 ENDIF652 !653 END SUBROUTINE654 655 656 527 SUBROUTINE tra_asm_inc( kt ) 657 528 !!---------------------------------------------------------------------- … … 664 535 !! ** Action : 665 536 !!---------------------------------------------------------------------- 666 INTEGER, INTENT(IN) :: kt! Current time step667 ! 668 INTEGER :: ji,jj,jk669 INTEGER :: it537 INTEGER, INTENT(IN) :: kt ! Current time step 538 ! 539 INTEGER :: ji, jj, jk 540 INTEGER :: it 670 541 REAL(wp) :: zincwgt ! IAU weight for current time step 671 542 REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 672 543 !!---------------------------------------------------------------------- 673 544 ! 674 545 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 675 546 ! used to prevent the applied increments taking the temperature below the local freezing point 676 677 DO jk=1, jpkm1 678 fzptnz (:,:,jk) = tfreez( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) ) 679 ENDDO 680 681 IF ( ln_asmiau ) THEN 682 683 !-------------------------------------------------------------------- 684 ! Incremental Analysis Updating 685 !-------------------------------------------------------------------- 686 547 DO jk = 1, jpkm1 548 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) ) 549 END DO 550 ! 551 ! !-------------------------------------- 552 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 553 ! !-------------------------------------- 554 ! 687 555 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 688 556 ! 689 557 it = kt - nit000 + 1 690 558 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 691 559 ! 692 560 IF(lwp) THEN 693 561 WRITE(numout,*) 694 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', & 695 & kt,' with IAU weight = ', wgtiau(it) 562 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 696 563 WRITE(numout,*) '~~~~~~~~~~~~' 697 564 ENDIF 698 565 ! 699 566 ! Update the tracer tendencies 700 567 DO jk = 1, jpkm1 … … 719 586 ENDIF 720 587 END DO 721 722 ENDIF 723 588 ! 589 ENDIF 590 ! 724 591 IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work 725 592 DEALLOCATE( t_bkginc ) 726 593 DEALLOCATE( s_bkginc ) 727 594 ENDIF 728 729 730 ELSEIF ( ln_asmdin ) THEN 731 732 !-------------------------------------------------------------------- 733 ! Direct Initialization 734 !-------------------------------------------------------------------- 735 595 ! !-------------------------------------- 596 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 597 ! !-------------------------------------- 598 ! 736 599 IF ( kt == nitdin_r ) THEN 737 600 ! 738 601 neuler = 0 ! Force Euler forward step 739 602 ! 740 603 ! Initialize the now fields with the background + increment 741 604 IF (ln_temnofreeze) THEN 742 605 ! Do not apply negative increments if the temperature will fall below freezing 743 WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 744 & tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 606 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 745 607 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 746 608 END WHERE … … 751 613 ! Do not apply negative increments if the salinity will fall below a specified 752 614 ! minimum value salfixmin 753 WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 754 & tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 615 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 755 616 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 756 617 END WHERE … … 759 620 ENDIF 760 621 761 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields 762 763 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities 764 765 IF( ln_zps .AND. .NOT. lk_c1d ) & 766 & CALL zps_hde( nit000, jpts, tsb, & ! Partial steps: before horizontal derivative 767 & gtsu, gtsv, rhd, & ! of T, S, rd at the bottom ocean level 768 & gru , grv ) 769 770 #if defined key_zdfkpp 771 CALL eos( tsn, rhd, fsdept_n(:,:,:) ) ! Compute rhd 772 #endif 622 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields 623 624 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities 625 !!gm fabien 626 ! CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities 627 !!gm 628 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 773 635 774 636 DEALLOCATE( t_bkginc ) … … 780 642 ENDIF 781 643 ! Perhaps the following call should be in step 782 IF 644 IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment 783 645 ! 784 646 END SUBROUTINE tra_asm_inc … … 801 663 REAL(wp) :: zincwgt ! IAU weight for current time step 802 664 !!---------------------------------------------------------------------- 803 804 IF ( ln_asmiau ) THEN 805 806 !-------------------------------------------------------------------- 807 ! Incremental Analysis Updating 808 !-------------------------------------------------------------------- 809 665 ! 666 ! !-------------------------------------------- 667 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 668 ! !-------------------------------------------- 669 ! 810 670 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 811 671 ! 812 672 it = kt - nit000 + 1 813 673 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 814 674 ! 815 675 IF(lwp) THEN 816 676 WRITE(numout,*) 817 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', & 818 & kt,' with IAU weight = ', wgtiau(it) 677 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 819 678 WRITE(numout,*) '~~~~~~~~~~~~' 820 679 ENDIF 821 680 ! 822 681 ! Update the dynamic tendencies 823 682 DO jk = 1, jpkm1 … … 825 684 va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt 826 685 END DO 827 686 ! 828 687 IF ( kt == nitiaufin_r ) THEN 829 688 DEALLOCATE( u_bkginc ) 830 689 DEALLOCATE( v_bkginc ) 831 690 ENDIF 832 833 ENDIF 834 835 ELSEIF ( ln_asmdin ) THEN 836 837 !-------------------------------------------------------------------- 838 ! Direct Initialization 839 !-------------------------------------------------------------------- 840 691 ! 692 ENDIF 693 ! !----------------------------------------- 694 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 695 ! !----------------------------------------- 696 ! 841 697 IF ( kt == nitdin_r ) THEN 842 698 ! 843 699 neuler = 0 ! Force Euler forward step 844 700 ! 845 701 ! Initialize the now fields with the background + increment 846 702 un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:) 847 703 vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 848 704 ! 849 705 ub(:,:,:) = un(:,:,:) ! Update before fields 850 706 vb(:,:,:) = vn(:,:,:) 851 707 ! 852 708 DEALLOCATE( u_bkg ) 853 709 DEALLOCATE( v_bkg ) … … 877 733 REAL(wp) :: zincwgt ! IAU weight for current time step 878 734 !!---------------------------------------------------------------------- 879 880 IF ( ln_asmiau ) THEN 881 882 !-------------------------------------------------------------------- 883 ! Incremental Analysis Updating 884 !-------------------------------------------------------------------- 885 735 ! 736 ! !----------------------------------------- 737 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 738 ! !----------------------------------------- 739 ! 886 740 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 887 741 ! 888 742 it = kt - nit000 + 1 889 743 zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step 890 744 ! 891 745 IF(lwp) THEN 892 746 WRITE(numout,*) … … 895 749 WRITE(numout,*) '~~~~~~~~~~~~' 896 750 ENDIF 897 751 ! 898 752 ! Save the tendency associated with the IAU weighted SSH increment 899 753 ! (applied in dynspg.*) … … 904 758 DEALLOCATE( ssh_bkginc ) 905 759 ENDIF 906 907 ENDIF 908 909 ELSEIF ( ln_asmdin ) THEN 910 911 !-------------------------------------------------------------------- 912 ! Direct Initialization 913 !-------------------------------------------------------------------- 914 760 ! 761 ENDIF 762 ! !----------------------------------------- 763 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 764 ! !----------------------------------------- 765 ! 915 766 IF ( kt == nitdin_r ) THEN 916 917 neuler = 0 ! Force Euler forward step 918 919 ! Initialize the now fields the background + increment 920 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 921 922 ! Update before fields 923 sshb(:,:) = sshn(:,:) 924 925 IF( lk_vvl ) THEN 926 DO jk = 1, jpk 927 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 928 END DO 929 ENDIF 930 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 ! 931 776 DEALLOCATE( ssh_bkg ) 932 777 DEALLOCATE( ssh_bkginc ) 933 778 ! 934 779 ENDIF 935 780 ! … … 950 795 !! 951 796 !!---------------------------------------------------------------------- 952 IMPLICIT NONE 953 ! 954 INTEGER, INTENT(in) :: kt ! Current time step 797 INTEGER, INTENT(in) :: kt ! Current time step 955 798 INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation 956 799 ! … … 962 805 #endif 963 806 !!---------------------------------------------------------------------- 964 965 IF ( ln_asmiau ) THEN 966 967 !-------------------------------------------------------------------- 968 ! Incremental Analysis Updating 969 !-------------------------------------------------------------------- 970 807 ! 808 ! !----------------------------------------- 809 IF ( ln_asmiau ) THEN ! Incremental Analysis Updating 810 ! !----------------------------------------- 811 ! 971 812 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 972 813 ! 973 814 it = kt - nit000 + 1 974 815 zincwgt = wgtiau(it) ! IAU weight for the current time step 975 816 ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 976 817 ! 977 818 IF(lwp) THEN 978 819 WRITE(numout,*) 979 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 980 & 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) 981 821 WRITE(numout,*) '~~~~~~~~~~~~' 982 822 ENDIF 983 823 ! 984 824 ! Sea-ice : LIM-3 case (to add) 985 825 ! 986 826 #if defined key_lim2 987 827 ! Sea-ice : LIM-2 case … … 1021 861 1022 862 #if defined key_cice && defined key_asminc 1023 ! Sea-ice : CICE case. Zero ice increment tendency into CICE 1024 ndaice_da(:,:) = 0.0_wp 1025 #endif 1026 1027 ENDIF 1028 1029 ELSEIF ( ln_asmdin ) THEN 1030 1031 !-------------------------------------------------------------------- 1032 ! Direct Initialization 1033 !-------------------------------------------------------------------- 1034 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 ! 1035 871 IF ( kt == nitdin_r ) THEN 1036 872 ! 1037 873 neuler = 0 ! Force Euler forward step 1038 874 ! 1039 875 ! Sea-ice : LIM-3 case (to add) 1040 876 ! 1041 877 #if defined key_lim2 1042 878 ! Sea-ice : LIM-2 case. … … 1054 890 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 1055 891 ELSEWHERE 1056 zhicifinc(:,:) = 0. 0_wp892 zhicifinc(:,:) = 0._wp 1057 893 END WHERE 1058 894 ! … … 1063 899 ! seaice salinity balancing (to add) 1064 900 #endif 1065 901 ! 1066 902 #if defined key_cice && defined key_asminc 1067 903 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 1068 904 ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 1069 905 #endif 1070 IF ( .NOT. PRESENT(kindic) ) THEN1071 DEALLOCATE( seaice_bkginc )1072 END IF1073 906 IF ( .NOT. PRESENT(kindic) ) THEN 907 DEALLOCATE( seaice_bkginc ) 908 END IF 909 ! 1074 910 ELSE 1075 911 ! 1076 912 #if defined key_cice && defined key_asminc 1077 ! Sea-ice : CICE case. Zero ice increment tendency into CICE1078 ndaice_da(:,:) = 0.0_wp 1079 #endif 1080 913 ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE 914 915 #endif 916 ! 1081 917 ENDIF 1082 918 … … 1155 991 ! 1156 992 !#endif 1157 993 ! 1158 994 ENDIF 1159 995 ! 1160 996 END SUBROUTINE seaice_asm_inc 1161 997
Note: See TracChangeset
for help on using the changeset viewer.
