- Timestamp:
- 2016-10-21T17:38:13+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90
r6140 r7068 41 41 PUBLIC trd_mxl_trc 42 42 PUBLIC trd_mxl_trc_alloc 43 PUBLIC trd_mxl_bio44 43 PUBLIC trd_mxl_trc_init 45 44 PUBLIC trd_mxl_trc_zint 46 PUBLIC trd_mxl_bio_zint47 45 48 46 CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file 49 47 INTEGER :: nmoymltrd 50 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndextrd1 51 INTEGER, DIMENSION(jptra) :: nidtrd, nh_t 48 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndextrd1, nidtrd, nh_t 52 49 INTEGER :: ndimtrd1 53 50 INTEGER, SAVE :: ionce, icount 54 #if defined key_pisces_reduced55 INTEGER :: nidtrdbio, nh_tb56 INTEGER, SAVE :: ioncebio, icountbio57 INTEGER, SAVE :: nmoymltrdbio58 #endif59 51 LOGICAL :: llwarn = .TRUE. ! this should always be .TRUE. 60 52 LOGICAL :: lldebug = .TRUE. 61 53 62 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: ztmltrd2 ! 63 #if defined key_pisces_reduced64 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ztmltrdbio2 ! only needed for mean diagnostics in trd_mxl_bio()65 #endif66 55 67 56 !! * Substitutions … … 79 68 !!---------------------------------------------------------------------- 80 69 ALLOCATE( ztmltrd2(jpi,jpj,jpltrd_trc,jptra) , & 81 #if defined key_pisces_reduced 82 & ztmltrdbio2(jpi,jpj,jpdiabio) , & 83 #endif 84 & ndextrd1(jpi*jpj) , STAT=trd_mxl_trc_alloc) 70 & ndextrd1(jpi*jpj), nidtrd(jptra), nh_t(jptra), STAT=trd_mxl_trc_alloc) 85 71 ! 86 72 IF( lk_mpp ) CALL mpp_sum ( trd_mxl_trc_alloc ) … … 131 117 SELECT CASE ( nn_ctls_trc ) ! choice of the control surface 132 118 CASE ( -2 ) ; STOP 'trdmxl_trc : not ready ' ! -> isopycnal surface (see ???) 133 #if defined key_pisces || defined key_pisces_reduced134 119 CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion 135 #endif136 120 CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl) 137 121 CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file … … 207 191 ! 208 192 END SUBROUTINE trd_mxl_trc_zint 209 210 211 SUBROUTINE trd_mxl_bio_zint( ptrc_trdmxl, ktrd )212 !!----------------------------------------------------------------------213 !! *** ROUTINE trd_mxl_bio_zint ***214 !!215 !! ** Purpose : Compute the vertical average of the 3D fields given as arguments216 !! to the subroutine. This vertical average is performed from ocean217 !! surface down to a chosen control surface.218 !!219 !! ** Method/usage :220 !! The control surface can be either a mixed layer depth (time varying)221 !! or a fixed surface (jk level or bowl).222 !! Choose control surface with nctls in namelist NAMTRD :223 !! nctls_trc = 0 : use mixed layer with density criterion224 !! nctls_trc = 1 : read index from file 'ctlsurf_idx'225 !! nctls_trc > 1 : use fixed level surface jk = nctls_trc226 !! Note: in the remainder of the routine, the volume between the227 !! surface and the control surface is called "mixed-layer"228 !!----------------------------------------------------------------------229 !!230 INTEGER , INTENT(in) :: ktrd ! bio trend index231 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptrc_trdmxl ! passive trc trend232 #if defined key_pisces_reduced233 !234 INTEGER :: ji, jj, jk, isum235 REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk236 !!----------------------------------------------------------------------237 238 CALL wrk_alloc( jpi, jpj, zvlmsk )239 240 ! I. Definition of control surface and integration weights241 ! --------------------------------------------------------242 ! ==> only once per time step <==243 244 IF( icountbio == 1 ) THEN245 !246 tmltrd_bio(:,:,:) = 0.e0 ! <<< reset trend arrays to zero247 ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer248 SELECT CASE ( nn_ctls_trc ) ! choice of the control surface249 CASE ( -2 ) ; STOP 'trdmxl_trc : not ready ' ! -> isopycnal surface (see ???)250 CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion251 CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl)252 CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file253 CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )254 nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level255 END SELECT256 257 ! ... Compute ndextrd1 and ndimtrd1 only once258 IF( ioncebio == 1 ) THEN259 !260 ! Check of validity : nmld_trc(ji,jj) <= jpktrd_trc261 isum = 0262 zvlmsk(:,:) = 0.e0263 264 IF( jpktrd_trc < jpk ) THEN265 DO jj = 1, jpj266 DO ji = 1, jpi267 IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN268 zvlmsk(ji,jj) = tmask(ji,jj,1)269 ELSE270 isum = isum + 1271 zvlmsk(ji,jj) = 0.272 END IF273 END DO274 END DO275 END IF276 277 ! Index of ocean points (2D only)278 IF( isum > 0 ) THEN279 WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum280 CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )281 ELSE282 CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )283 END IF284 285 ioncebio = 0 ! no more pass here286 !287 END IF ! ( ioncebio == 1 )288 289 ! ... Weights for vertical averaging290 wkx_trc(:,:,:) = 0.e0291 DO jk = 1, jpktrd_trc ! initialize wkx_trc with vertical scale factor in mixed-layer292 DO jj = 1,jpj293 DO ji = 1,jpi294 IF( jk - nmld_trc(ji,jj) < 0. ) wkx_trc(ji,jj,jk) = e3t_n(ji,jj,jk) * tmask(ji,jj,jk)295 END DO296 END DO297 END DO298 299 rmld_trc(:,:) = 0.300 DO jk = 1, jpktrd_trc ! compute mixed-layer depth : rmld_trc301 rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)302 END DO303 304 DO jk = 1, jpktrd_trc ! compute integration weights305 wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )306 END DO307 308 icountbio = 0 ! <<< flag = off : control surface & integr. weights309 ! ! computed only once per time step310 END IF ! ( icountbio == 1 )311 312 ! II. Vertical integration of trends in the mixed-layer313 ! -----------------------------------------------------314 315 316 DO jk = 1, jpktrd_trc317 tmltrd_bio(:,:,ktrd) = tmltrd_bio(:,:,ktrd) + ptrc_trdmxl(:,:,jk) * wkx_trc(:,:,jk)318 END DO319 320 CALL wrk_dealloc( jpi, jpj, zvlmsk )321 #endif322 !323 END SUBROUTINE trd_mxl_bio_zint324 193 325 194 … … 877 746 878 747 879 SUBROUTINE trd_mxl_bio( kt )880 !!----------------------------------------------------------------------881 !! *** ROUTINE trd_mld ***882 !!883 !! ** Purpose : Compute and cumulate the mixed layer biological trends over an analysis884 !! period, and write NetCDF outputs.885 !!886 !! ** Method/usage :887 !! The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant888 !! logical namelist variable) :889 !! 1) to explain the difference between initial and final890 !! mixed-layer T & S (where initial and final relate to the891 !! current analysis window, defined by ntrd in the namelist)892 !! 2) to explain the difference between the current and previous893 !! TIME-AVERAGED mixed-layer T & S (where time-averaging is894 !! performed over each analysis window).895 !!896 !! ** Consistency check :897 !! If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt898 !! entrainment) should be zero, at machine accuracy. Note that in the case899 !! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO900 !! over the first two analysis windows (except if restart).901 !! N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8902 !! for checking residuals.903 !! On a NEC-SX5 computer, this typically leads to:904 !! O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.905 !! O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.906 !!907 !! ** Action :908 !! At each time step, mixed-layer averaged trends are stored in the909 !! tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).910 !! This array is known when trd_mld is called, at the end of the stp subroutine,911 !! except for the purely vertical K_z diffusion term, which is embedded in the912 !! lateral diffusion trend.913 !!914 !! In I), this K_z term is diagnosed and stored, thus its contribution is removed915 !! from the lateral diffusion trend.916 !! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative917 !! arrays are updated.918 !! In III), called only once per analysis window, we compute the total trends,919 !! along with the residuals and the Asselin correction terms.920 !! In IV), the appropriate trends are written in the trends NetCDF file.921 !!922 !! References :923 !! - Vialard & al.924 !! - See NEMO documentation (in preparation)925 !!----------------------------------------------------------------------926 INTEGER, INTENT( in ) :: kt ! ocean time-step index927 #if defined key_pisces_reduced928 INTEGER :: jl, it, itmod929 LOGICAL :: llwarn = .TRUE., lldebug = .TRUE.930 REAL(wp) :: zfn, zfn2931 !!----------------------------------------------------------------------932 ! ... Warnings933 IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )934 935 ! ======================================================================936 ! II. Cumulate the trends over the analysis window937 ! ======================================================================938 939 ztmltrdbio2(:,:,:) = 0.e0 ! <<< reset arrays to zero940 941 ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window942 ! ------------------------------------------------------------------------943 IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)944 !945 tmltrd_csum_ub_bio (:,:,:) = 0.e0946 !947 END IF948 949 ! II.4 Cumulated trends over the analysis period950 ! ----------------------------------------------951 !952 ! [ 1rst analysis window ] [ 2nd analysis window ]953 !954 !955 ! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps956 ! ntrd 2*ntrd etc.957 ! 1 2 3 4 =5 e.g. =10958 !959 IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN960 !961 nmoymltrdbio = nmoymltrdbio + 1962 963 ! ... Trends associated with the time mean of the ML passive tracers964 tmltrd_sum_bio (:,:,:) = tmltrd_sum_bio (:,:,:) + tmltrd_bio (:,:,:)965 tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)966 !967 END IF968 969 ! ======================================================================970 ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)971 ! ======================================================================972 973 ! Convert to appropriate physical units974 tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * rn_ucf_trc975 976 MODULO_NTRD : IF( MOD( kt, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd977 !978 zfn = float(nmoymltrdbio) ; zfn2 = zfn * zfn979 980 ! III.1 Prepare fields for output ("instantaneous" diagnostics)981 ! -------------------------------------------------------------982 983 #if defined key_diainstant984 STOP 'tmltrd_bio : key_diainstant was never checked within trdmxl. Comment this to proceed.'985 #endif986 ! III.2 Prepare fields for output ("mean" diagnostics)987 ! ----------------------------------------------------988 989 ztmltrdbio2(:,:,:) = tmltrd_csum_ub_bio(:,:,:) + tmltrd_csum_ln_bio(:,:,:)990 991 !-- Lateral boundary conditions992 IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration993 ! ES_B27_CD_WARN : lbc inutile GYRE, cf. + haut994 DO jn = 1, jpdiabio995 CALL lbc_lnk( ztmltrdbio2(:,:,jn), 'T', 1. )996 ENDDO997 ENDIF998 999 IF( lldebug ) THEN1000 !1001 WRITE(numout,*) 'trd_mxl_bio : write trends in the Mixed Layer for debugging process:'1002 WRITE(numout,*) '~~~~~~~~~~~ '1003 WRITE(numout,*) 'TRC kt = ', kt, 'nmoymltrdbio = ', nmoymltrdbio1004 WRITE(numout,*)1005 1006 DO jl = 1, jpdiabio1007 IF( ln_trdmxl_trc_instant ) THEN1008 WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &1009 & ' SUM tmltrd_bio : ', SUM2D(tmltrd_bio(:,:,jl))1010 ELSE1011 WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &1012 & ' SUM ztmltrdbio2 : ', SUM2D(ztmltrdbio2(:,:,jl))1013 endif1014 END DO1015 1016 97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)1017 98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)1018 99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))1019 WRITE(numout,*)1020 !1021 ENDIF1022 1023 ! III.3 Time evolution array swap1024 ! -------------------------------1025 1026 ! For passive tracer mean diagnostics1027 tmltrd_csum_ub_bio (:,:,:) = zfn * tmltrd_sum_bio(:,:,:) - tmltrd_csum_ln_bio(:,:,:)1028 1029 ! III.4 Convert to appropriate physical units1030 ! -------------------------------------------1031 ztmltrdbio2 (:,:,:) = ztmltrdbio2 (:,:,:) * rn_ucf_trc/zfn21032 1033 END IF MODULO_NTRD1034 1035 ! ======================================================================1036 ! IV. Write trends in the NetCDF file1037 ! ======================================================================1038 1039 ! IV.1 Code for IOIPSL/NetCDF output1040 ! ----------------------------------1041 1042 ! define time axis1043 itmod = kt - nittrc000 + 11044 it = kt1045 1046 IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN1047 WRITE(numout,*) ' '1048 WRITE(numout,*) 'trd_mxl_bio : write ML bio trends in the NetCDF file :'1049 WRITE(numout,*) '~~~~~~~~~~~ '1050 WRITE(numout,*) ' ', TRIM(clhstnam), ' at kt = ', kt1051 WRITE(numout,*) ' N.B. nmoymltrdbio = ', nmoymltrdbio1052 WRITE(numout,*) ' '1053 END IF1054 1055 1056 ! 2. Start writing data1057 ! ---------------------1058 1059 NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags1060 !1061 DO jl = 1, jpdiabio1062 CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) , &1063 & it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )1064 END DO1065 1066 1067 IF( kt == nitend ) CALL histclo( nidtrdbio )1068 1069 ELSE ! <<< write the trends for passive tracer mean diagnostics1070 1071 DO jl = 1, jpdiabio1072 CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) , &1073 & it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )1074 END DO1075 1076 IF( kt == nitend ) CALL histclo( nidtrdbio )1077 !1078 END IF NETCDF_OUTPUT1079 1080 ! Compute the control surface (for next time step) : flag = on1081 icountbio = 11082 1083 1084 1085 IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN1086 !1087 ! III.5 Reset cumulative arrays to zero1088 ! -------------------------------------1089 nmoymltrdbio = 01090 tmltrd_csum_ln_bio (:,:,:) = 0.e01091 tmltrd_sum_bio (:,:,:) = 0.e01092 END IF1093 1094 ! ======================================================================1095 ! Write restart file1096 ! ======================================================================1097 1098 ! restart write is done in trd_mxl_trc_write which is called by trd_mxl_bio (Marina)1099 !1100 #endif1101 END SUBROUTINE trd_mxl_bio1102 1103 1104 748 REAL FUNCTION sum2d( ztab ) 1105 749 !!---------------------------------------------------------------------- … … 1191 835 tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0 1192 836 1193 #if defined key_pisces_reduced1194 nmoymltrdbio = 01195 tmltrd_sum_bio (:,:,:) = 0.e0 ; tmltrd_csum_ln_bio (:,:,:) = 0.e01196 DO jl = 1, jp_pisces_trd1197 ctrd_bio(jl,1) = ctrbil(jl) ! long name1198 ctrd_bio(jl,2) = ctrbio(jl) ! short name1199 ENDDO1200 #endif1201 1202 837 IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN 1203 838 CALL trd_mxl_trc_rst_read … … 1208 843 tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean 1209 844 tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0 1210 #if defined key_pisces_reduced1211 tmltrd_csum_ub_bio (:,:,:) = 0.e01212 #endif1213 845 1214 846 ENDIF … … 1216 848 icount = 1 ; ionce = 1 ! open specifier 1217 849 1218 #if defined key_pisces_reduced1219 icountbio = 1 ; ioncebio = 1 ! open specifier1220 #endif1221 850 1222 851 ! I.3 Read control surface from file ctlsurf_idx … … 1308 937 END DO 1309 938 1310 #if defined key_pisces_reduced1311 !-- Create a NetCDF file and enter the define mode1312 CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' )1313 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &1314 & 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom, snc4chunks=snc4set )1315 #endif1316 1317 939 !-- Define physical units 1318 940 IF( rn_ucf_trc == 1. ) THEN … … 1354 976 END DO 1355 977 1356 #if defined key_pisces_reduced1357 DO jl = 1, jp_pisces_trd1358 CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1)) , &1359 & cltrcu, jpi, jpj, nh_tb, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean1360 END DO ! if zsto=rdt above1361 #endif1362 1363 978 !-- Leave IOIPSL/NetCDF define mode 1364 979 DO jn = 1, jptra … … 1366 981 END DO 1367 982 1368 #if defined key_pisces_reduced1369 !-- Leave IOIPSL/NetCDF define mode1370 CALL histend( nidtrdbio, snc4set )1371 1372 983 IF(lwp) WRITE(numout,*) 1373 IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'1374 #endif1375 984 1376 985 END SUBROUTINE trd_mxl_trc_init … … 1385 994 WRITE(*,*) 'trd_mxl_trc: You should not have seen this print! error?', kt 1386 995 END SUBROUTINE trd_mxl_trc 1387 SUBROUTINE trd_mxl_bio( kt )1388 INTEGER, INTENT( in) :: kt1389 WRITE(*,*) 'trd_mxl_bio: You should not have seen this print! error?', kt1390 END SUBROUTINE trd_mxl_bio1391 996 SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn ) 1392 997 INTEGER , INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
Note: See TracChangeset
for help on using the changeset viewer.
