Changeset 14072 for NEMO/trunk/src/OCE/ASM/asminc.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/ASM/asminc.F90 (modified) (42 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ASM/asminc.F90
r13982 r14072 9 9 !! ! 2007-04 (A. Weaver) Merge with OPAVAR/NEMOVAR 10 10 !! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 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 13 !! ! 2014-09 (D. Lea) Local calc_date removed use routine from OBS … … 32 32 USE zpshde ! Partial step : Horizontal Derivative 33 33 USE asmpar ! Parameters for the assmilation interface 34 USE asmbkg ! 34 USE asmbkg ! 35 35 USE c1d ! 1D initialization 36 36 USE sbc_oce ! Surface boundary condition variables. … … 46 46 IMPLICIT NONE 47 47 PRIVATE 48 48 49 49 PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights 50 50 PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments … … 73 73 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u- & v- velocity components 74 74 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S 75 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 75 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 76 76 REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step 77 77 #if defined key_asminc … … 81 81 INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term 82 82 INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization 83 INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval 83 INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval 84 84 INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval 85 ! 85 ! 86 86 INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting 87 ! !: = 1 Linear hat-like, centred in middle of IAU interval 87 ! !: = 1 Linear hat-like, centred in middle of IAU interval 88 88 REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix) 89 89 … … 107 107 !!---------------------------------------------------------------------- 108 108 !! *** ROUTINE asm_inc_init *** 109 !! 109 !! 110 110 !! ** Purpose : Initialize the assimilation increment and IAU weights. 111 111 !! 112 112 !! ** Method : Initialize the assimilation increment and IAU weights. 113 113 !! 114 !! ** Action : 114 !! ** Action : 115 115 !!---------------------------------------------------------------------- 116 116 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices … … 264 264 ! 265 265 ! !--------------------------------------------------------- 266 IF( niaufn == 0 ) THEN ! Constant IAU forcing 266 IF( niaufn == 0 ) THEN ! Constant IAU forcing 267 267 ! !--------------------------------------------------------- 268 268 DO jt = 1, iiauper … … 270 270 END DO 271 271 ! !--------------------------------------------------------- 272 ELSEIF ( niaufn == 1 ) THEN ! Linear hat-like, centred in middle of IAU interval 272 ELSEIF ( niaufn == 1 ) THEN ! Linear hat-like, centred in middle of IAU interval 273 273 ! !--------------------------------------------------------- 274 274 ! Compute the normalization factor 275 275 znorm = 0._wp 276 276 IF( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval 277 imid = iiauper / 2 277 imid = iiauper / 2 278 278 DO jt = 1, imid 279 279 znorm = znorm + REAL( jt ) … … 281 281 znorm = 2.0 * znorm 282 282 ELSE ! Odd number of time steps in IAU interval 283 imid = ( iiauper + 1 ) / 2 283 imid = ( iiauper + 1 ) / 2 284 284 DO jt = 1, imid - 1 285 285 znorm = znorm + REAL( jt ) … … 308 308 DO jt = 1, icycper 309 309 ztotwgt = ztotwgt + wgtiau(jt) 310 WRITE(numout,*) ' ', jt, ' ', wgtiau(jt) 311 END DO 310 WRITE(numout,*) ' ', jt, ' ', wgtiau(jt) 311 END DO 312 312 WRITE(numout,*) ' ===================================' 313 313 WRITE(numout,*) ' Time-integrated weight = ', ztotwgt 314 314 WRITE(numout,*) ' ===================================' 315 315 ENDIF 316 316 317 317 ENDIF 318 318 … … 339 339 CALL iom_open( c_asminc, inum ) 340 340 ! 341 CALL iom_get( inum, 'time' , zdate_inc ) 341 CALL iom_get( inum, 'time' , zdate_inc ) 342 342 CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb ) 343 343 CALL iom_get( inum, 'z_inc_datef', z_inc_datef ) … … 346 346 ! 347 347 IF(lwp) THEN 348 WRITE(numout,*) 348 WRITE(numout,*) 349 349 WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef 350 350 WRITE(numout,*) '~~~~~~~~~~~~' … … 360 360 & ' not agree with Direct Initialization time' ) 361 361 362 IF ( ln_trainc ) THEN 362 IF ( ln_trainc ) THEN 363 363 CALL iom_get( inum, jpdom_auto, 'bckint', t_bkginc, 1 ) 364 364 CALL iom_get( inum, jpdom_auto, 'bckins', s_bkginc, 1 ) … … 372 372 ENDIF 373 373 374 IF ( ln_dyninc ) THEN 375 CALL iom_get( inum, jpdom_auto, 'bckinu', u_bkginc, 1 ) 376 CALL iom_get( inum, jpdom_auto, 'bckinv', v_bkginc, 1 ) 374 IF ( ln_dyninc ) THEN 375 CALL iom_get( inum, jpdom_auto, 'bckinu', u_bkginc, 1 ) 376 CALL iom_get( inum, jpdom_auto, 'bckinv', v_bkginc, 1 ) 377 377 ! Apply the masks 378 378 u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:) … … 383 383 WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0 384 384 ENDIF 385 385 386 386 IF ( ln_sshinc ) THEN 387 387 CALL iom_get( inum, jpdom_auto, 'bckineta', ssh_bkginc, 1 ) … … 409 409 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN ! Apply divergence damping filter 410 410 ! !-------------------------------------- 411 ALLOCATE( zhdiv(jpi,jpj) ) 411 ALLOCATE( zhdiv(jpi,jpj) ) 412 412 ! 413 413 DO jt = 1, nn_divdmp … … 428 428 & + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 429 429 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) & 430 & + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 430 & + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 431 431 END_2D 432 432 END DO … … 434 434 END DO 435 435 ! 436 DEALLOCATE( zhdiv ) 436 DEALLOCATE( zhdiv ) 437 437 ! 438 438 ENDIF … … 455 455 CALL iom_open( c_asmdin, inum ) 456 456 ! 457 CALL iom_get( inum, 'rdastp', zdate_bkg ) 457 CALL iom_get( inum, 'rdastp', zdate_bkg ) 458 458 ! 459 459 IF(lwp) THEN 460 WRITE(numout,*) 460 WRITE(numout,*) 461 461 WRITE(numout,*) ' ==>>> Assimilation background state valid at : ', zdate_bkg 462 462 WRITE(numout,*) … … 467 467 & ' not agree with Direct Initialization time' ) 468 468 ! 469 IF ( ln_trainc ) THEN 469 IF ( ln_trainc ) THEN 470 470 CALL iom_get( inum, jpdom_auto, 'tn', t_bkg ) 471 471 CALL iom_get( inum, jpdom_auto, 'sn', s_bkg ) … … 474 474 ENDIF 475 475 ! 476 IF ( ln_dyninc ) THEN 476 IF ( ln_dyninc ) THEN 477 477 CALL iom_get( inum, jpdom_auto, 'un', u_bkg, cd_type = 'U', psgn = 1._wp ) 478 478 CALL iom_get( inum, jpdom_auto, 'vn', v_bkg, cd_type = 'V', psgn = 1._wp ) … … 502 502 ! 503 503 END SUBROUTINE asm_inc_init 504 505 504 505 506 506 SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs ) 507 507 !!---------------------------------------------------------------------- 508 508 !! *** ROUTINE tra_asm_inc *** 509 !! 509 !! 510 510 !! ** Purpose : Apply the tracer (T and S) assimilation increments 511 511 !! 512 512 !! ** Method : Direct initialization or Incremental Analysis Updating 513 513 !! 514 !! ** Action : 514 !! ** Action : 515 515 !!---------------------------------------------------------------------- 516 516 INTEGER , INTENT(in ) :: kt ! Current time step … … 524 524 !!---------------------------------------------------------------------- 525 525 ! 526 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 527 ! used to prevent the applied increments taking the temperature below the local freezing point 526 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 527 ! used to prevent the applied increments taking the temperature below the local freezing point 528 528 IF( ln_temnofreeze ) THEN 529 529 DO jk = 1, jpkm1 … … 587 587 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 588 588 ! !-------------------------------------- 589 ! 589 ! 590 590 IF ( kt == nitdin_r ) THEN 591 591 ! … … 647 647 ! 648 648 ENDIF 649 ! 649 ! 650 650 ENDIF 651 651 ! Perhaps the following call should be in step … … 658 658 !!---------------------------------------------------------------------- 659 659 !! *** ROUTINE dyn_asm_inc *** 660 !! 660 !! 661 661 !! ** Purpose : Apply the dynamics (u and v) assimilation increments. 662 662 !! 663 663 !! ** Method : Direct initialization or Incremental Analysis Updating. 664 664 !! 665 !! ** Action : 665 !! ** Action : 666 666 !!---------------------------------------------------------------------- 667 667 INTEGER , INTENT( in ) :: kt ! ocean time-step index … … 684 684 ! 685 685 IF(lwp) THEN 686 WRITE(numout,*) 686 WRITE(numout,*) 687 687 WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 688 688 WRITE(numout,*) '~~~~~~~~~~~~' … … 704 704 ELSEIF ( ln_asmdin ) THEN ! Direct Initialization 705 705 ! !----------------------------------------- 706 ! 706 ! 707 707 IF ( kt == nitdin_r ) THEN 708 708 ! … … 711 711 ! Initialize the now fields with the background + increment 712 712 puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:) 713 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 713 pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 714 714 ! 715 715 puu(:,:,:,Kbb) = puu(:,:,:,Kmm) ! Update before fields … … 730 730 !!---------------------------------------------------------------------- 731 731 !! *** ROUTINE ssh_asm_inc *** 732 !! 732 !! 733 733 !! ** Purpose : Apply the sea surface height assimilation increment. 734 734 !! 735 735 !! ** Method : Direct initialization or Incremental Analysis Updating. 736 736 !! 737 !! ** Action : 737 !! ** Action : 738 738 !!---------------------------------------------------------------------- 739 739 INTEGER, INTENT(IN) :: kt ! Current time step … … 755 755 ! 756 756 IF(lwp) THEN 757 WRITE(numout,*) 757 WRITE(numout,*) 758 758 WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', & 759 759 & kt,' with IAU weight = ', wgtiau(it) … … 807 807 !! *** ROUTINE ssh_asm_div *** 808 808 !! 809 !! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence 809 !! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence 810 810 !! across all the water column 811 811 !! … … 823 823 REAL(wp), DIMENSION(:,:) , POINTER :: ztim ! local array 824 824 !!---------------------------------------------------------------------- 825 ! 825 ! 826 826 #if defined key_asminc 827 827 CALL ssh_asm_inc( kt, Kbb, Kmm ) !== (calculate increments) 828 828 ! 829 IF( ln_linssh ) THEN 829 IF( ln_linssh ) THEN 830 830 phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1) 831 ELSE 831 ELSE 832 832 ALLOCATE( ztim(jpi,jpj) ) 833 833 ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) ) 834 DO jk = 1, jpkm1 835 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 834 DO jk = 1, jpkm1 835 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 836 836 END DO 837 837 ! … … 846 846 !!---------------------------------------------------------------------- 847 847 !! *** ROUTINE seaice_asm_inc *** 848 !! 848 !! 849 849 !! ** Purpose : Apply the sea ice assimilation increment. 850 850 !! 851 851 !! ** Method : Direct initialization or Incremental Analysis Updating. 852 852 !! 853 !! ** Action : 853 !! ** Action : 854 854 !! 855 855 !!---------------------------------------------------------------------- … … 873 873 ! 874 874 it = kt - nit000 + 1 875 zincwgt = wgtiau(it) ! IAU weight for the current time step 875 zincwgt = wgtiau(it) ! IAU weight for the current time step 876 876 ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments) 877 877 ! … … 997 997 !#if defined defined key_si3 || defined key_cice 998 998 ! 999 ! IF (ln_seaicebal ) THEN 999 ! IF (ln_seaicebal ) THEN 1000 1000 ! !! balancing salinity increments 1001 1001 ! !! simple case from limflx.F90 (doesn't include a mass flux) … … 1009 1009 ! 1010 1010 ! DO jj = 1, jpj 1011 ! DO ji = 1, jpi 1011 ! DO ji = 1, jpi 1012 1012 ! ! calculate change in ice and snow mass per unit area 1013 1013 ! ! positive values imply adding salt to the ocean (results from ice formation) … … 1020 1020 ! 1021 1021 ! ! prevent small mld 1022 ! ! less than 10m can cause salinity instability 1022 ! ! less than 10m can cause salinity instability 1023 1023 ! IF (mld < 10) mld=10 1024 1024 ! 1025 ! ! set to bottom of a level 1025 ! ! set to bottom of a level 1026 1026 ! DO jk = jpk-1, 2, -1 1027 1027 ! IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN … … 1032 1032 ! 1033 1033 ! ! avoid applying salinity balancing in shallow water or on land 1034 ! ! 1034 ! ! 1035 1035 ! 1036 1036 ! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) … … 1043 1043 ! 1044 1044 ! ! put increments in for levels in the mixed layer 1045 ! ! but prevent salinity below a threshold value 1046 ! 1047 ! DO jk = 1, jkmax 1048 ! 1049 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1045 ! ! but prevent salinity below a threshold value 1046 ! 1047 ! DO jk = 1, jkmax 1048 ! 1049 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1050 1050 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 1051 1051 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn … … 1058 1058 ! ! 1059 1059 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean 1060 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1061 ! !! 1062 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1060 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1061 ! !! 1062 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1063 1063 ! !! ! E-P (kg m-2 s-2) 1064 1064 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) … … 1073 1073 ! 1074 1074 END SUBROUTINE seaice_asm_inc 1075 1075 1076 1076 !!====================================================================== 1077 1077 END MODULE asminc
Note: See TracChangeset
for help on using the changeset viewer.
