Changeset 5099 for LMDZ6/branches/Amaury_dev/libf/phylmd/cosp
- Timestamp:
- Jul 22, 2024, 9:29:09 PM (4 months ago)
- Location:
- LMDZ6/branches/Amaury_dev/libf/phylmd/cosp
- Files:
-
- 36 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/MISR_simulator.F
r5095 r5099 1 ! 1 2 2 ! Copyright (c) 2009, Roger Marchand, version 1.2 3 3 ! All rights reserved. 4 4 ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ 5 5 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MISR_simulator/MISR_simulator.f $ 6 ! 6 7 7 ! Redistribution and use in source and binary forms, with or without modification, are permitted 8 8 ! provided that the following conditions are met: 9 ! 9 10 10 ! * Redistributions of source code must retain the above copyright notice, this list of 11 11 ! conditions and the following disclaimer. … … 15 15 ! * Neither the name of the University of Washington nor the names of its contributors may be used 16 16 ! to endorse or promote products derived from this software without specific prior written permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 19 19 ! BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT … … 22 22 ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 23 23 ! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 24 !25 24 26 25 SUBROUTINE MISR_simulator( … … 119 118 120 119 tauchk = -1.*log(0.9999999) 121 122 ! 120 123 121 ! For each GCM cell or horizontal model grid point ... 124 ! 122 125 123 do j=1,npoints 126 124 127 !128 125 ! estimate distribution of Model layer tops 129 ! 126 130 127 dist_model_layertops(j,:)=0 131 128 … … 154 151 & dist_model_layertops(j,iMISR_ztop)+1 155 152 enddo 156 157 158 ! 153 159 154 ! compute total cloud optical depth for each column 160 ! 155 161 156 do ibox=1,ncol 162 157 … … 280 275 enddo ! loop of subcolumns 281 276 enddo ! loop of gridpoints 282 283 284 ! 277 285 278 ! Modify MISR CTH for satellite spatial / pattern matcher effects 286 ! 279 287 280 ! Code in this region added by roj 5/2006 to account 288 281 ! for spatial effect of the MISR pattern matcher. … … 291 284 ! a lower CTH, THEN misr will tend to but place the 292 285 ! odd column at the same height as it neighbors. 293 ! 286 294 287 ! This setup assumes the columns represent a about a 1 to 4 km scale 295 288 ! it will need to be modified significantly, otherwise … … 337 330 endif 338 331 339 !340 332 ! DETERMINE CLOUD TYPE FREQUENCIES 341 ! 333 342 334 ! Now that ztop and tau have been determined, 343 335 ! determine amount of each cloud type … … 407 399 408 400 else 409 410 ! 401 411 402 ! determine index for MISR bin set 412 !413 403 414 404 iMISR_ztop=2 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/array_lib.F90
r5095 r5099 4 4 ! Compiled/Modified: 5 5 ! 07/01/06 John Haynes (haynes@atmos.colostate.edu) 6 ! 6 7 7 ! infind (function) 8 8 ! lin_interpolate (function) … … 19 19 use m_mrgrnk 20 20 implicit none 21 ! 21 22 22 ! Purpose: 23 23 ! Finds the index of an array that is closest to a value, plus the 24 24 ! difference between the value found and the value specified 25 ! 25 26 26 ! Inputs: 27 27 ! [list] an array of sequential values … … 29 29 ! Optional input: 30 30 ! [sort] set to 1 if [list] is in unknown/non-sequential order 31 ! 31 32 32 ! Returns: 33 33 ! index of [list] that is closest to [val] 34 ! 34 35 35 ! Optional output: 36 36 ! [dist] set to variable containing [list([result])] - [val] 37 ! 37 38 38 ! Requires: 39 39 ! mrgrnk library 40 ! 40 41 41 ! Created: 42 42 ! 10/16/03 John Haynes (haynes@atmos.colostate.edu) … … 101 101 use m_mrgrnk 102 102 implicit none 103 ! 103 104 104 ! Purpose: 105 105 ! linearly interpolate a set of y2 values given a set of y1,x1,x2 106 ! 106 107 107 ! Inputs: 108 108 ! [yarr] an array of y1 values … … 110 110 ! [xxarr] an array of x2 values 111 111 ! [tol] maximum distance for a match 112 ! 112 113 113 ! Output: 114 114 ! [yyarr] interpolated array of y2 values 115 ! 115 116 116 ! Requires: 117 117 ! mrgrnk library 118 ! 118 119 119 ! Created: 120 120 ! 06/07/06 John Haynes (haynes@atmos.colostate.edu) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/atmos_lib.F90
r5095 r5099 4 4 ! Compiled/Modified: 5 5 ! 07/01/06 John Haynes (haynes@atmos.colostate.edu) 6 ! 6 7 7 ! mcclatchey (subroutine) 8 8 … … 17 17 subroutine mcclatchey(stype,hgt,prs,tk,rh) 18 18 implicit none 19 ! 19 20 20 ! Purpose: 21 21 ! returns a standard atmospheric profile 22 ! 22 23 23 ! Input: 24 24 ! [stype] type of profile to return … … 26 26 ! 2 = mid-latitude winter 27 27 ! 3 = tropical 28 ! 28 29 29 ! Outputs: 30 30 ! [hgt] height (m) … … 32 32 ! [tk] temperature (K) 33 33 ! [rh] relative humidity (%) 34 ! 34 35 35 ! Created: 36 36 ! 06/01/2006 John Haynes (haynes@atmos.colostate.edu) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/calc_Re.F90
r5095 r5099 7 7 ! Purpose: 8 8 ! Calculates Effective Radius (1/2 distribution 3rd moment / 2nd moment). 9 ! 9 10 10 ! For some distribution types, the total number concentration (per kg), Np 11 11 ! may be optionally specified. Should be set to zero, otherwise. 12 ! 12 13 13 ! Roj Marchand July 2010 14 14 15 15 16 16 ! Inputs: 17 ! 17 18 18 ! [Q] hydrometeor mixing ratio (g/kg) ! not needed for some distribution types 19 19 ! [Np] Optional Total number concentration (per kg). 0 = use defaults (p1, p2, p3) 20 ! 20 21 21 ! [rho_a] ambient air density (kg m^-3) 22 ! 22 23 23 ! Distribution parameters as per quickbeam documentation. 24 24 ! [dtype] distribution type … … 28 28 ! [bmp] b params for mass 29 29 ! [p1],[p2],[p3] distribution parameters 30 ! 31 ! 30 31 32 32 ! Outputs: 33 33 ! [Re] Effective radius, 1/2 the 3rd moment/2nd moment (um) 34 ! 34 35 35 ! Created: 36 36 ! July 2010 Roj Marchand 37 ! 38 39 37 40 38 ! ----- INPUTS ----- 41 39 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/congvec.h
r2428 r5099 5 5 ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ 6 6 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/icarus-scops-4.1-bsd/congvec.f $ 7 ! 7 8 8 ! Redistribution and use in source and binary forms, with or without 9 9 ! modification, are permitted provided that the 10 10 ! following conditions are met: 11 ! 11 12 12 ! * Redistributions of source code must retain the above 13 13 ! copyright notice, this list of conditions and the following … … 21 21 ! derived from this software without specific prior written 22 22 ! permission. 23 ! 23 24 24 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 25 25 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT … … 33 33 ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 34 34 ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 35 ! 35 36 36 ! *****************************COPYRIGHT******************************* 37 37 ! *****************************COPYRIGHT******************************* -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/cosp_defs.h
r2431 r5099 3 3 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_defs.h $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! to endorse or promote products derived from this software without specific prior written 16 16 ! permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/dsd.F90
r5095 r5099 1 ! 1 2 2 ! $Id$ 3 ! 3 4 4 subroutine dsd(Q,Re_,Np,D,N,nsizes,dtype,rho_a,tk, & 5 5 dmin,dmax,apm,bpm,rho_c,p1,p2,p3) … … 10 10 ! Purpose: 11 11 ! Create a discrete drop size distribution 12 ! 12 13 13 ! Starting with Quickbeam V3, this routine now allows input of 14 14 ! both effective radius (Re) and total number concentration (Nt) 15 15 ! Roj Marchand July 2010 16 ! 16 17 17 ! The version in Quickbeam v.104 was modified to allow Re but not Nt 18 18 ! This is a significantly modified form for the version 19 ! 19 20 20 ! Originally Part of QuickBeam v1.03 by John Haynes 21 21 ! http://reef.atmos.colostate.edu/haynes/radarsim 22 ! 22 23 23 ! Inputs: 24 ! 24 25 25 ! [Q] hydrometeor mixing ratio (g/kg) 26 26 ! [Re] Optional Effective Radius (microns). 0 = use defaults (p1, p2, p3) 27 ! 27 28 28 ! [D] array of discrete drop sizes (um) where we desire to know the number concentraiton n(D). 29 29 ! [nsizes] number of elements of [D] 30 ! 30 31 31 ! [dtype] distribution type 32 32 ! [rho_a] ambient air density (kg m^-3) … … 36 36 ! [rho_c] alternate constant density (kg m^-3) 37 37 ! [p1],[p2],[p3] distribution parameters 38 ! 38 39 39 ! Input/Output: 40 40 ! [apm] a parameter for mass (kg m^[-bpm]) 41 41 ! [bmp] b params for mass 42 ! 42 43 43 ! Outputs: 44 44 ! [N] discrete concentrations (cm^-3 um^-1) 45 45 ! or, for monodisperse, a constant (1/cm^3) 46 ! 46 47 47 ! Requires: 48 48 ! function infind 49 ! 49 50 50 ! Created: 51 51 ! 11/28/05 John Haynes (haynes@atmos.colostate.edu) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/format_input.F90
r5095 r5099 4 4 ! Compiled/Modified: 5 5 ! 08/28/2006 John Haynes (haynes@atmos.colostate.edu) 6 ! 6 7 7 ! irreg_to_grid (subroutine) 8 8 ! order_data (subroutine) … … 23 23 ! Linearly interpolate sounding-level data to the hydrometeor-level 24 24 ! resolution 25 ! 25 26 26 ! Inputs: 27 27 ! [hgt_matrix] hydrometeor-level heights … … 30 30 ! [env_p_matrix] sounding-level pressures 31 31 ! [env_rh_matrix] sounding-level relative humidities 32 ! 32 33 33 ! Outputs: 34 34 ! [t_matrix] hydrometeor-level temperatures 35 35 ! [p_matrix] hydrometeor-level pressures 36 36 ! [rh_matrix] hydrometeor-level relative humidities 37 ! 37 38 38 ! Created: 39 39 ! 08/28/2006 John Haynes (haynes@atmos.colostate.edu) … … 73 73 ! Ensure that input data is in top-down order/bottom-up order, 74 74 ! for space-based/surface based radars, respectively 75 ! 75 76 76 ! Inputs: 77 77 ! [hgt_matrix] heights … … 81 81 ! [rh_matrix] relative humidities 82 82 ! [sfc_radar] 1=surface radar, 0=spaceborne 83 ! 83 84 84 ! Outputs: 85 85 ! [hgt_matrix],[hm_matrix],[p_matrix,[t_matrix],[rh_matrix] in proper 86 86 ! order for input to the radar simulator routine 87 87 ! [hgt_reversed] T=heights were reordered,F=heights were not reordered 88 ! 88 89 89 ! Note: 90 90 ! The order for all profiles is assumed to the same as the first profile. 91 ! 91 92 92 ! Created: 93 93 ! 08/28/2006 John Haynes (haynes@atmos.colostate.edu) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/gases.F90
r5095 r5099 6 6 ! Purpose: 7 7 ! Compute 2-way gaseous attenuation through a volume in microwave 8 ! 8 9 9 ! Inputs: 10 10 ! [PRES_mb] pressure (mb) (hPa) … … 12 12 ! [RH] relative humidity (%) 13 13 ! [f] frequency (GHz), < 300 GHz 14 ! 14 15 15 ! Returns: 16 16 ! 2-way gaseous attenuation (dB/km) 17 ! 17 18 18 ! Reference: 19 19 ! Uses method of Liebe (1985) 20 ! 20 21 21 ! Created: 22 22 ! 12/09/05 John Haynes (haynes@atmos.colostate.edu) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/icarus.F
r5095 r5099 41 41 ! Corporation. 42 42 ! All rights reserved. 43 ! 43 44 44 ! Redistribution and use in source and binary forms, with or without 45 45 ! modification, are permitted provided that the 46 46 ! following conditions are met: 47 ! 47 48 48 ! * Redistributions of source code must retain the above 49 49 ! copyright notice, this list of conditions and the following … … 58 58 ! derived from this software without specific prior written 59 59 ! permission. 60 ! 60 61 61 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 62 62 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT … … 70 70 ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 71 71 ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 72 ! 72 73 73 ! *****************************COPYRIGHT******************************* 74 74 ! *****************************COPYRIGHT******************************* … … 153 153 ! with interpolated temperature equal to the radiance 154 154 ! determined cloud-top temperature 155 ! 155 156 156 ! 1 = find the *lowest* altitude (highest pressure) level 157 157 ! with interpolated temperature equal to the radiance 158 158 ! determined cloud-top temperature 159 ! 159 160 160 ! 2 = find the *highest* altitude (lowest pressure) level 161 161 ! with interpolated temperature equal to the radiance 162 162 ! determined cloud-top temperature 163 ! 163 164 164 ! ONLY APPLICABLE IF top_height EQUALS 1 or 3 165 165 ! ! 166 166 ! 1 = old setting: matches all versions of 167 167 ! ISCCP simulator with versions numbers 3.5.1 and lower 168 ! 168 169 169 ! 2 = default setting: for version numbers 4.0 and higher 170 ! 170 171 171 ! The following input variables are used only if top_height = 1 or top_height = 3 172 ! 172 173 173 REAL skt(npoints) ! skin Temperature (K) 174 174 REAL emsfc_lw ! 10.5 micron emissivity of surface (fraction) … … 222 222 223 223 REAL boxptop(npoints,ncol) ! cloud top pressure (mb) in each column 224 225 226 ! 224 227 225 ! ------ 228 226 ! Working variables added when program updated to mimic Mark Webb's PV-Wave code … … 456 454 enddo 457 455 458 !459 456 ! ---------------------------------------------------! 460 457 461 462 !463 458 ! ---------------------------------------------------! 464 459 ! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and … … 502 497 enddo 503 498 endif 504 ! 499 505 500 ! ---------------------------------------------------! 506 501 507 508 509 !510 502 ! ---------------------------------------------------! 511 503 ! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES 512 504 ! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE 513 ! 505 514 506 ! again this is only done if top_height = 1 or 3 515 ! 507 516 508 ! fluxtop is the 10.5 micron radiance at the top of the 517 509 ! atmosphere … … 525 517 526 518 !---------------------------------------------------------------------- 527 ! 519 528 520 ! DO CLEAR SKY RADIANCE CALCULATION FIRST 529 ! 521 530 522 !compute water vapor continuum emissivity 531 523 !this treatment follows Schwarkzopf and Ramasamy … … 640 632 enddo 641 633 endif 642 643 644 ! 634 645 635 ! END OF CLEAR SKY CALCULATION 646 ! 636 647 637 !---------------------------------------------------------------- 648 638 … … 795 785 !clouds which have partial emissivity and need the 796 786 !adjustment performed in this section 797 ! 787 798 788 !If it turns out that the cloud brightness temperature 799 789 !is greater than 260K, then the liquid cloud conversion 800 790 !factor of 2.56 is used. 801 ! 791 802 792 !Note that this is discussed on pages 85-87 of 803 793 !the ISCCP D level documentation (Rossow et al. 1996) … … 915 905 ! ---------------------------------------------------! 916 906 917 !918 907 ! ---------------------------------------------------! 919 908 ! DETERMINE CLOUD TOP PRESSURE 920 ! 909 921 910 ! again the 2 methods differ according to whether 922 911 ! or not you use the physical cloud top pressure (top_height = 2) 923 912 ! or the radiatively determined cloud top pressure (top_height = 1 or 3) 924 !925 913 926 914 !compute cloud top pressure … … 1004 992 1005 993 30 continue 1006 1007 ! 1008 ! 994 995 1009 996 ! ---------------------------------------------------! 1010 997 1011 1012 !1013 998 ! ---------------------------------------------------! 1014 999 ! DETERMINE ISCCP CLOUD TYPE FREQUENCIES 1015 ! 1000 1016 1001 ! Now that ptop and tau have been determined, 1017 1002 ! determine amount of each of the 49 ISCCP cloud 1018 1003 ! types 1019 ! 1004 1020 1005 ! Also compute grid box mean cloud top pressure and 1021 1006 ! optical thickness. The mean cloud top pressure and … … 1027 1012 ! then converting the average albedo back to a mean 1028 1013 ! optical thickness. 1029 !1030 1014 1031 1015 !compute isccp frequencies … … 1187 1171 end if 1188 1172 enddo ! j 1189 ! 1173 1190 1174 ! ---------------------------------------------------! 1191 1175 1192 1176 ! ---------------------------------------------------! 1193 1177 ! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM 1194 ! 1178 1195 1179 if (debugcol.ne.0) then 1196 ! 1180 1197 1181 do j=1,npoints,debugcol 1198 1182 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/isccp_cloud_types.F
r2432 r5099 41 41 ! (c) British Crown Copyright 2009, the Met Office. 42 42 ! All rights reserved. 43 ! 43 44 44 ! Redistribution and use in source and binary forms, with or without 45 45 ! modification, are permitted provided that the 46 46 ! following conditions are met: 47 ! 47 48 48 ! * Redistributions of source code must retain the above 49 49 ! copyright notice, this list of conditions and the following … … 57 57 ! derived from this software without specific prior written 58 58 ! permission. 59 ! 59 60 60 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 61 61 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT … … 69 69 ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 70 70 ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 71 ! 71 72 72 ! *****************************COPYRIGHT******************************* 73 73 ! *****************************COPYRIGHT******************************* … … 161 161 ! with interpolated temperature equal to the radiance 162 162 ! determined cloud-top temperature 163 ! 163 164 164 ! 1 = find the *lowest* altitude (highest pressure) level 165 165 ! with interpolated temperature equal to the radiance 166 166 ! determined cloud-top temperature 167 ! 167 168 168 ! 2 = find the *highest* altitude (lowest pressure) level 169 169 ! with interpolated temperature equal to the radiance 170 170 ! determined cloud-top temperature 171 ! 171 172 172 ! ONLY APPLICABLE IF top_height EQUALS 1 or 3 173 ! 173 174 174 ! 1 = old setting: matches all versions of 175 175 ! ISCCP simulator with versions numbers 3.5.1 and lower 176 ! 176 177 177 ! 2 = default setting: for version numbers 4.0 and higher 178 ! 178 179 179 ! The following input variables are used only if top_height = 1 or top_height = 3 180 ! 180 181 181 REAL skt(npoints) ! skin Temperature (K) 182 182 REAL emsfc_lw ! 10.5 micron emissivity of surface (fraction) … … 230 230 231 231 REAL boxptop(npoints,ncol) ! cloud top pressure (mb) in each column 232 233 234 ! 232 235 233 ! ------ 236 234 ! Working variables added when program updated to mimic Mark Webb's PV-Wave code -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/lidar_simulator.F90
r5095 r5099 3 3 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/actsim/lidar_simulator.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! contributors may be used to endorse or promote products derived from this software without 16 16 ! specific prior written permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 31 31 , ls_radliq, ls_radice, cv_radliq, cv_radice & 32 32 , ice_type, pmol, pnorm, pnorm_perp_tot,tautot, refl ) 33 ! 33 34 34 !--------------------------------------------------------------------------------- 35 35 ! Purpose: To compute lidar signal from model-simulated profiles of cloud water 36 36 ! and cloud fraction in each sub-column of each model gridbox. 37 ! 37 38 38 ! References: 39 39 ! Chepfer H., S. Bony, D. Winker, M. Chiriaco, J.-L. Dufresne, G. Seze (2008), 40 40 ! Use of CALIPSO lidar observations to evaluate the cloudiness simulated by a 41 41 ! climate model, Geophys. Res. Lett., 35, L15704, doi:10.1029/2008GL034207. 42 ! 42 43 43 ! Previous references: 44 44 ! Chiriaco et al, MWR, 2006; Chepfer et al., MWR, 2007 45 ! 45 46 46 ! Contacts: Helene Chepfer (chepfer@lmd.polytechnique.fr), Sandrine Bony (bony@lmd.jussieu.fr) 47 ! 47 48 48 ! May 2007: ActSim code of M. Chiriaco and H. Chepfer rewritten by S. Bony 49 ! 49 50 50 ! May 2008, H. Chepfer: 51 51 ! - Units of pressure inputs: Pa 52 52 ! - Non Spherical particles : LS Ice NS coefficients, CONV Ice NS coefficients 53 53 ! - New input: ice_type (0=ice-spheres ; 1=ice-non-spherical) 54 ! 54 55 55 ! June 2008, A. Bodas-Salcedo: 56 56 ! - Ported to Fortran 90 and optimisation changes 57 ! 57 58 58 ! August 2008, J-L Dufresne: 59 59 ! - Optimisation changes (sum instructions suppressed) 60 ! 60 61 61 ! October 2008, S. Bony, H. Chepfer and J-L. Dufresne : 62 62 ! - Interface with COSP v2.0: … … 65 65 ! depolarisation diagnostic removed 66 66 ! parasol (polder) reflectances (for 5 different solar zenith angles) added 67 ! 67 68 68 ! December 2008, S. Bony, H. Chepfer and J-L. Dufresne : 69 69 ! - Modification of the integration of the lidar equation. 70 70 ! - change the cloud detection threshold 71 ! 71 72 72 ! April 2008, A. Bodas-Salcedo: 73 73 ! - Bug fix in computation of pmol and pnorm of upper layer 74 ! 74 75 75 ! April 2008, J-L. Dufresne 76 76 ! - Bug fix in computation of pmol and pnorm, thanks to Masaki Satoh: a factor 2 77 77 ! was missing. This affects the ATB values but not the cloud fraction. 78 ! 78 79 79 ! January 2013, G. Cesana and H. Chepfer: 80 80 ! - Add the perpendicular component of the backscattered signal (pnorm_perp_tot) in the arguments … … 83 83 ! Reference: Cesana G. and H. Chepfer (2013): Evaluation of the cloud water phase 84 84 ! in a climate model using CALIPSO-GOCCP, J. Geophys. Res., doi: 10.1002/jgrd.50376 85 ! 85 86 86 !--------------------------------------------------------------------------------- 87 ! 87 88 88 ! Inputs: 89 89 ! npoints : number of horizontal points … … 106 106 ! ice_type : ice particle shape hypothesis (ice_type=0 for spheres, ice_type=1 107 107 ! for non spherical particles) 108 ! 108 109 109 ! Outputs: 110 110 ! pmol : molecular attenuated backscatter lidar signal power (m^-1.sr^-1) … … 113 113 ! tautot: optical thickess integrated from top to level z 114 114 ! refl : parasol(polder) reflectance 115 ! 115 116 116 ! Version 1.0 (June 2007) 117 117 ! Version 1.1 (May 2008) … … 228 228 ! as a function of the ATB for ice or liquid cloud particles derived from CALIPSO-GOCCP 229 229 ! observations at 120m vertical grid (Cesana and Chepfer, JGR, 2013). 230 ! 230 231 231 ! Relationship between ATBice and ATBperp,ice for ice particles 232 232 ! ATBperp,ice = Alpha*ATBice … … 414 414 415 415 ! Total signal (molecular + particules): 416 ! 417 ! 416 417 418 418 ! For performance reason on vector computers, the 2 following lines should not be used 419 419 ! and should be replace by the later one. … … 426 426 tautot(:,:) = tautot(:,:) + tau_part(:,:,i) 427 427 ENDDO ! i 428 ! 428 429 429 ! Upper layer 430 430 pnorm(:,nlev) = betatot(:,nlev) / (2.*tautot(:,nlev)) & … … 585 585 !--------------------------------------------------------- 586 586 ! Parasol/Polder module 587 ! 587 588 588 ! Purpose : Compute reflectance for one particular viewing direction 589 589 ! and 5 solar zenith angles (calculation valid only over ocean) … … 611 611 612 612 END SUBROUTINE lidar_simulator 613 ! 613 614 614 !--------------------------------------------------------------------------------- 615 ! 615 616 616 SUBROUTINE parasol(npoints,nrefl,undef & 617 617 ,tautot_S_liq,tautot_S_ice & … … 621 621 ! of cloud water and cloud fraction in each sub-column of each model 622 622 ! gridbox. 623 ! 624 ! 623 624 625 625 ! December 2008, S. Bony, H. Chepfer and J-L. Dufresne : 626 626 ! - optimization for vectorization 627 ! 627 628 628 ! Version 2.0 (October 2008) 629 629 ! Version 2.1 (December 2008) … … 642 642 ! outputs 643 643 REAL refl(npoints,nrefl) ! Parasol reflectances 644 ! 644 645 645 ! Local variables 646 646 REAL tautot_S(npoints) ! cloud optical thickness, from TOA to surface … … 702 702 rlumA_mod=0 703 703 rlumB_mod=0 704 ! 704 705 705 pi = ACOS(-1.0) 706 706 r_norm(:)=1./ COS(pi/180.*tetas(:)) 707 ! 707 708 708 tautot_S_liq(:)=MAX(tautot_S_liq(:),tau(1)) 709 709 tautot_S_ice(:)=MAX(tautot_S_ice(:),tau(1)) 710 710 tautot_S(:) = tautot_S_ice(:) + tautot_S_liq(:) 711 ! 711 712 712 ! relative fraction of the opt. thick due to liquid or ice clouds 713 713 WHERE (tautot_S(:) .GT. 0.) … … 719 719 END WHERE 720 720 tautot_S(:)=MIN(tautot_S(:),tau(nbtau)) 721 ! 721 722 722 ! Linear interpolation : 723 723 … … 730 730 bB(:,ny) = rlumB(:,ny) - aB(:,ny)*tau(ny) 731 731 ENDDO 732 ! 732 733 733 DO it=1,ntetas 734 734 DO ny=1,nbtau-1 … … 739 739 END DO 740 740 END DO 741 ! 741 742 742 DO it=1,ntetas 743 743 refl(:,it) = frac_taucol_liq(:) * rlumA_mod(:,it) & -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/m_mrgrnk.F90
r5095 r5099 21 21 ! __________________________________________________________ 22 22 Real (kind=kdp) :: XVALA, XVALB 23 ! 23 24 24 Integer, Dimension (SIZE(IRNGT)) :: JWRKT 25 25 Integer :: LMTNA, LMTNC, IRNG1, IRNG2 26 26 Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB 27 ! 27 28 28 NVAL = Min (SIZE(XDONT), SIZE(IRNGT)) 29 29 Select Case (NVAL) … … 36 36 Continue 37 37 End Select 38 ! 38 39 39 ! Fill-in the index array, creating ordered couples 40 ! 40 41 41 Do IIND = 2, NVAL, 2 42 42 If (XDONT(IIND-1) <= XDONT(IIND)) Then … … 51 51 IRNGT (NVAL) = NVAL 52 52 End If 53 ! 53 54 54 ! We will now have ordered subsets A - B - A - B - ... 55 55 ! and merge A and B couples into C - C - ... 56 ! 56 57 57 LMTNA = 2 58 58 LMTNC = 4 59 ! 59 60 60 ! First iteration. The length of the ordered subsets goes from 2 to 4 61 ! 61 62 62 Do 63 63 If (NVAL <= 2) Exit 64 ! 64 65 65 ! Loop on merges of A and B into C 66 ! 66 67 67 Do IWRKD = 0, NVAL - 1, 4 68 68 If ((IWRKD+4) > NVAL) Then 69 69 If ((IWRKD+2) >= NVAL) Exit 70 ! 70 71 71 ! 1 2 3 72 ! 72 73 73 If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit 74 ! 74 75 75 ! 1 3 2 76 ! 76 77 77 If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then 78 78 IRNG2 = IRNGT (IWRKD+2) 79 79 IRNGT (IWRKD+2) = IRNGT (IWRKD+3) 80 80 IRNGT (IWRKD+3) = IRNG2 81 ! 81 82 82 ! 3 1 2 83 ! 83 84 84 Else 85 85 IRNG1 = IRNGT (IWRKD+1) … … 90 90 Exit 91 91 End If 92 ! 92 93 93 ! 1 2 3 4 94 ! 94 95 95 If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle 96 ! 96 97 97 ! 1 3 x x 98 ! 98 99 99 If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then 100 100 IRNG2 = IRNGT (IWRKD+2) … … 108 108 IRNGT (IWRKD+4) = IRNG2 109 109 End If 110 ! 110 111 111 ! 3 x x x 112 ! 112 113 113 Else 114 114 IRNG1 = IRNGT (IWRKD+1) … … 133 133 End If 134 134 End Do 135 ! 135 136 136 ! The Cs become As and Bs 137 ! 137 138 138 LMTNA = 4 139 139 Exit 140 140 End Do 141 ! 141 142 142 ! Iteration loop. Each time, the length of the ordered subsets 143 143 ! is doubled. 144 ! 144 145 145 Do 146 146 If (LMTNA >= NVAL) Exit 147 147 IWRKF = 0 148 148 LMTNC = 2 * LMTNC 149 ! 149 150 150 ! Loop on merges of A and B into C 151 ! 151 152 152 Do 153 153 IWRK = IWRKF … … 161 161 IINDA = 1 162 162 IINDB = JINDA + 1 163 ! 163 164 164 ! Shortcut for the case when the max of A is smaller 165 165 ! than the min of B. This line may be activated when the 166 166 ! initial set is already close to sorted. 167 ! 167 168 168 ! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE 169 ! 169 170 170 ! One steps in the C subset, that we build in the final rank array 171 ! 171 172 172 ! Make a copy of the rank array for the merge iteration 173 ! 173 174 174 JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA) 175 ! 175 176 176 XVALA = XDONT (JWRKT(IINDA)) 177 177 XVALB = XDONT (IRNGT(IINDB)) 178 ! 178 179 179 Do 180 180 IWRK = IWRK + 1 181 ! 181 182 182 ! We still have unprocessed values in both A and B 183 ! 183 184 184 If (XVALA > XVALB) Then 185 185 IRNGT (IWRK) = IRNGT (IINDB) … … 197 197 XVALA = XDONT (JWRKT(IINDA)) 198 198 End If 199 ! 199 200 200 End Do 201 201 End Do 202 ! 202 203 203 ! The Cs become As and Bs 204 ! 204 205 205 LMTNA = 2 * LMTNA 206 206 End Do 207 ! 207 208 208 Return 209 ! 209 210 210 End Subroutine D_mrgrnk 211 211 … … 221 221 ! __________________________________________________________ 222 222 Integer :: XVALA, XVALB 223 ! 223 224 224 Integer, Dimension (SIZE(IRNGT)) :: JWRKT 225 225 Integer :: LMTNA, LMTNC, IRNG1, IRNG2 226 226 Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB 227 ! 227 228 228 NVAL = Min (SIZE(XDONT), SIZE(IRNGT)) 229 229 Select Case (NVAL) … … 236 236 Continue 237 237 End Select 238 ! 238 239 239 ! Fill-in the index array, creating ordered couples 240 ! 240 241 241 Do IIND = 2, NVAL, 2 242 242 If (XDONT(IIND-1) <= XDONT(IIND)) Then … … 251 251 IRNGT (NVAL) = NVAL 252 252 End If 253 ! 253 254 254 ! We will now have ordered subsets A - B - A - B - ... 255 255 ! and merge A and B couples into C - C - ... 256 ! 256 257 257 LMTNA = 2 258 258 LMTNC = 4 259 ! 259 260 260 ! First iteration. The length of the ordered subsets goes from 2 to 4 261 ! 261 262 262 Do 263 263 If (NVAL <= 2) Exit 264 ! 264 265 265 ! Loop on merges of A and B into C 266 ! 266 267 267 Do IWRKD = 0, NVAL - 1, 4 268 268 If ((IWRKD+4) > NVAL) Then 269 269 If ((IWRKD+2) >= NVAL) Exit 270 ! 270 271 271 ! 1 2 3 272 ! 272 273 273 If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit 274 ! 274 275 275 ! 1 3 2 276 ! 276 277 277 If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then 278 278 IRNG2 = IRNGT (IWRKD+2) 279 279 IRNGT (IWRKD+2) = IRNGT (IWRKD+3) 280 280 IRNGT (IWRKD+3) = IRNG2 281 ! 281 282 282 ! 3 1 2 283 ! 283 284 284 Else 285 285 IRNG1 = IRNGT (IWRKD+1) … … 290 290 Exit 291 291 End If 292 ! 292 293 293 ! 1 2 3 4 294 ! 294 295 295 If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle 296 ! 296 297 297 ! 1 3 x x 298 ! 298 299 299 If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then 300 300 IRNG2 = IRNGT (IWRKD+2) … … 308 308 IRNGT (IWRKD+4) = IRNG2 309 309 End If 310 ! 310 311 311 ! 3 x x x 312 ! 312 313 313 Else 314 314 IRNG1 = IRNGT (IWRKD+1) … … 333 333 End If 334 334 End Do 335 ! 335 336 336 ! The Cs become As and Bs 337 ! 337 338 338 LMTNA = 4 339 339 Exit 340 340 End Do 341 ! 341 342 342 ! Iteration loop. Each time, the length of the ordered subsets 343 343 ! is doubled. 344 ! 344 345 345 Do 346 346 If (LMTNA >= NVAL) Exit 347 347 IWRKF = 0 348 348 LMTNC = 2 * LMTNC 349 ! 349 350 350 ! Loop on merges of A and B into C 351 ! 351 352 352 Do 353 353 IWRK = IWRKF … … 361 361 IINDA = 1 362 362 IINDB = JINDA + 1 363 ! 363 364 364 ! Shortcut for the case when the max of A is smaller 365 365 ! than the min of B. This line may be activated when the 366 366 ! initial set is already close to sorted. 367 ! 367 368 368 ! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE 369 ! 369 370 370 ! One steps in the C subset, that we build in the final rank array 371 ! 371 372 372 ! Make a copy of the rank array for the merge iteration 373 ! 373 374 374 JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA) 375 ! 375 376 376 XVALA = XDONT (JWRKT(IINDA)) 377 377 XVALB = XDONT (IRNGT(IINDB)) 378 ! 378 379 379 Do 380 380 IWRK = IWRK + 1 381 ! 381 382 382 ! We still have unprocessed values in both A and B 383 ! 383 384 384 If (XVALA > XVALB) Then 385 385 IRNGT (IWRK) = IRNGT (IINDB) … … 397 397 XVALA = XDONT (JWRKT(IINDA)) 398 398 End If 399 ! 399 400 400 End Do 401 401 End Do 402 ! 402 403 403 ! The Cs become As and Bs 404 ! 404 405 405 LMTNA = 2 * LMTNA 406 406 End Do 407 ! 407 408 408 Return 409 ! 409 410 410 End Subroutine I_mrgrnk 411 411 end module m_mrgrnk -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/math_lib.F90
r5095 r5099 4 4 ! Compiled/Modified: 5 5 ! 07/01/06 John Haynes (haynes@atmos.colostate.edu) 6 ! 6 7 7 ! gamma (function) 8 8 ! path_integral (function) … … 19 19 function gamma(x) 20 20 implicit none 21 ! 21 22 22 ! Purpose: 23 23 ! Returns the gamma function 24 ! 24 25 25 ! Input: 26 26 ! [x] value to compute gamma function of 27 ! 27 28 28 ! Returns: 29 29 ! gamma(x) 30 ! 30 31 31 ! Coded: 32 32 ! 02/02/06 John Haynes (haynes@atmos.colostate.edu) … … 100 100 use array_lib 101 101 implicit none 102 ! 102 103 103 ! Purpose: 104 104 ! evalues the integral (f ds) between f(index=i1) and f(index=i2) 105 105 ! using the AVINT procedure 106 ! 106 107 107 ! Inputs: 108 108 ! [f] functional values … … 110 110 ! [i1] index of lower limit 111 111 ! [i2] index of upper limit 112 ! 112 113 113 ! Returns: 114 114 ! result of path integral 115 ! 115 116 116 ! Notes: 117 117 ! [s] may be in forward or reverse numerical order 118 ! 118 119 119 ! Requires: 120 120 ! mrgrnk package 121 ! 121 122 122 ! Created: 123 123 ! 02/02/06 John Haynes (haynes@atmos.colostate.edu) … … 162 162 subroutine avint ( ftab, xtab, ntab, a_in, b_in, result ) 163 163 implicit none 164 ! 164 165 165 ! Purpose: 166 166 ! estimate the integral of unevenly spaced data 167 ! 167 168 168 ! Inputs: 169 169 ! [ftab] functional values … … 172 172 ! [a] lower limit of integration 173 173 ! [b] upper limit of integration 174 ! 174 175 175 ! Outputs: 176 176 ! [result] approximate value of integral 177 ! 177 178 178 ! Reference: 179 179 ! From SLATEC libraries, in public domain 180 ! 180 181 181 !*********************************************************************** 182 ! 182 183 183 ! AVINT estimates the integral of unevenly spaced data. 184 ! 184 185 185 ! Discussion: 186 ! 186 187 187 ! The method uses overlapping parabolas and smoothing. 188 ! 188 189 189 ! Modified: 190 ! 190 191 191 ! 30 October 2000 192 192 ! 4 January 2008, A. Bodas-Salcedo. Error control for XTAB taken out of 193 193 ! loop to allow vectorization. 194 ! 194 195 195 ! Reference: 196 ! 196 197 197 ! Philip Davis and Philip Rabinowitz, 198 198 ! Methods of Numerical Integration, 199 199 ! Blaisdell Publishing, 1967. 200 ! 200 201 201 ! P E Hennion, 202 202 ! Algorithm 77, … … 204 204 ! Communications of the Association for Computing Machinery, 205 205 ! Volume 5, page 96, 1962. 206 ! 206 207 207 ! Parameters: 208 ! 208 209 209 ! Input, real ( kind = 8 ) FTAB(NTAB), the function values, 210 210 ! FTAB(I) = F(XTAB(I)). 211 ! 211 212 212 ! Input, real ( kind = 8 ) XTAB(NTAB), the abscissas at which the 213 213 ! function values are given. The XTAB's must be distinct 214 214 ! and in ascending order. 215 ! 215 216 216 ! Input, integer NTAB, the number of entries in FTAB and 217 217 ! XTAB. NTAB must be at least 3. 218 ! 218 219 219 ! Input, real ( kind = 8 ) A, the lower limit of integration. A should 220 220 ! be, but need not be, near one endpoint of the interval 221 221 ! (X(1), X(NTAB)). 222 ! 222 223 223 ! Input, real ( kind = 8 ) B, the upper limit of integration. B should 224 224 ! be, but need not be, near one endpoint of the interval 225 225 ! (X(1), X(NTAB)). 226 ! 226 227 227 ! Output, real ( kind = 8 ) RESULT, the approximate value of the integral. 228 228 … … 293 293 return 294 294 end if 295 ! 295 296 296 ! If B < A, temporarily switch A and B, and store sign. 297 ! 297 298 298 if ( b < a ) then 299 299 syl = b … … 305 305 ind = 1 306 306 end if 307 ! 307 308 308 ! Bracket A and B between XTAB(ILO) and XTAB(IHI). 309 ! 309 310 310 ilo = 1 311 311 ihi = ntab … … 330 330 ihi = min ( ihi, ntab - 1 ) 331 331 ihi = max ( ilo, ihi - 1 ) 332 ! 332 333 333 ! Carry out approximate integration from XTAB(ILO) to XTAB(IHI). 334 ! 334 335 335 sum1 = 0.0D+00 336 336 … … 380 380 + cb * 0.5D+00 * ( b**2 - syl**2 ) & 381 381 + cc * ( b - syl ) 382 ! 382 383 383 ! Restore original values of A and B, reverse sign of integral 384 384 ! because of earlier switch. 385 ! 385 386 386 if ( ind /= 1 ) then 387 387 ind = 1 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp.F90
r5095 r5099 1 1 ! (c) British Crown Copyright 2008, the Met Office. 2 2 ! All rights reserved. 3 ! 3 4 4 ! Redistribution and use in source and binary forms, with or without modification, are permitted 5 5 ! provided that the following conditions are met: 6 ! 6 7 7 ! * Redistributions of source code must retain the above copyright notice, this list 8 8 ! of conditions and the following disclaimer. … … 13 13 ! to endorse or promote products derived from this software without specific prior written 14 14 ! permission. 15 ! 15 16 16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 17 17 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_constants.F90
r4785 r5099 1 1 ! (c) British Crown Copyright 2008, the Met Office. 2 2 ! All rights reserved. 3 ! 3 4 4 ! Redistribution and use in source and binary forms, with or without modification, are permitted 5 5 ! provided that the following conditions are met: 6 ! 6 7 7 ! * Redistributions of source code must retain the above copyright notice, this list 8 8 ! of conditions and the following disclaimer. … … 13 13 ! to endorse or promote products derived from this software without specific prior written 14 14 ! permission. 15 ! 15 16 16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 17 17 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 23 23 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 24 24 25 !26 25 ! History: 27 26 ! Jul 2007 - A. Bodas-Salcedo - Initial version … … 29 28 ! Oct 2008 - H. Chepfer - Added PARASOL_NREFL 30 29 ! Jun 2010 - R. Marchand - Modified to support quickbeam V3, added ifdef for hydrometeor definitions 31 ! 32 ! 33 ! 30 31 34 32 35 33 !!#INCLUDE "cosp_defs.h" … … 143 141 shape=(/2,MISR_N_CTH/)) 144 142 145 146 !147 143 ! The following code was modifed by Roj with implementation of quickbeam V3 148 144 ! (1) use ifdef to support more than one microphyscis scheme 149 145 ! (2) added constants microphysic_scheme_name, LOAD_scale_LUTs, and SAVE_scale_LUTs 150 !151 146 152 147 ! directory where LUTs will be stored … … 155 150 #ifdef MMF_V3_SINGLE_MOMENT 156 151 157 !158 152 ! Table hclass for quickbeam to support one-moment (bulk) microphysics scheme used by MMF V3.0 & V3.5 159 ! 160 161 ! 153 162 154 ! NOTE: if ANY value in this section of code is changed, the existing LUT 163 155 ! (i.e., the associated *.dat file) MUST be deleted so that a NEW 164 156 ! LUT will be created !!! 165 ! 157 166 158 character*120 :: RADAR_SIM_MICROPHYSICS_SCHEME_NAME = 'MMF_v3_single_moment' 167 159 … … 190 182 191 183 ! NOTES on HCLASS variables 192 ! 184 193 185 ! TYPE - Set to 194 186 ! 1 for modified gamma distribution, … … 217 209 ! P1, P2, P3 - are default distribution parameters that depend on the type 218 210 ! of distribution (see quickmbeam documentation for more information) 219 ! 211 220 212 ! Modified Gamma (must set P3 and one of P1 or P2) 221 213 ! P1 - Set to the total particle number concentration Nt /rho_a (kg-1 ), where … … 223 215 ! P2 - Set to the particle mean diameter D (micron). 224 216 ! P3 - Set to the distribution width nu. 225 ! 217 226 218 ! Exponetial (set one of) 227 219 ! P1 - Set to a constant intercept parameter N0 (m-4). 228 220 ! P2 - Set to a constant lambda (micron-1). 229 ! 221 230 222 ! Power Law 231 223 ! P1 - Set this to the value of a constant power law parameter br 232 ! 224 233 225 ! Monodisperse 234 226 ! P1 - Set to a constant diameter D0 (micron) = Re. 235 ! 227 236 228 ! Log-normal (must set P3 and one of P1 or P2) 237 229 ! P1 - Set to the total particle number concentration Nt /rho_a (kg-1 ) 238 230 ! P2 - Set to the geometric mean particle radius rg (micron). 239 231 ! P3 - Set to the natural logarithm of the geometric standard deviation. 240 !241 242 232 243 233 real,dimension(N_HYDRO) :: N_ax,N_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma_1,gamma_2,gamma_3,gamma_4 … … 265 255 #ifdef MMF_V3p5_TWO_MOMENT 266 256 267 !268 257 ! Table hclass for quickbeam to support two-moment "morrison" microphysics scheme used by V3.5 (SAM 6.8) 269 ! 258 270 259 ! This Number concentriation Np in [1/kg] MUST be input to COSP/radar simulator 271 ! 260 272 261 ! NOTE: Be sure to check that the ice-density (rho) set it this tables matches what you used 273 ! 274 275 ! 262 276 263 ! NOTE: if ANY value in this section of code is changed, the existing LUT 277 264 ! (i.e., the associated *.dat file) MUST be deleted so that a NEW 278 265 ! LUT will be created !!! 279 ! 266 280 267 character*120 :: RADAR_SIM_MICROPHYSICS_SCHEME_NAME = 'MMF_v3.5_two_moment' 281 268 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_isccp_simulator.F90
r3233 r5099 3 3 ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_isccp_simulator.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! to endorse or promote products derived from this software without specific prior written 16 16 ! permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_lidar.F90
r3233 r5099 3 3 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_lidar.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! to endorse or promote products derived from this software without specific prior written 16 16 ! permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 25 25 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 26 27 !28 27 ! History: 29 28 ! Jul 2007 - A. Bodas-Salcedo - Initial version … … 32 31 ! frac_out changed in sgx%frac_out) 33 32 ! Jun 2011 - G. Cesana - Added betaperp_tot argument 34 ! 35 ! 33 34 36 35 MODULE MOD_COSP_LIDAR 37 36 USE MOD_COSP_CONSTANTS -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_misr_simulator.F90
r3233 r5099 3 3 ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_misr_simulator.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! to endorse or promote products derived from this software without specific prior written 16 16 ! permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 25 25 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 26 27 !28 27 ! History: 29 28 ! Nov 2008 - A. Bodas-Salcedo - Initial version 30 ! 31 ! 29 32 30 33 31 MODULE MOD_COSP_MISR_SIMULATOR -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_modis_simulator.F90
r5095 r5099 4 4 ! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $ 5 5 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_modis_simulator.F90 $ 6 ! 6 7 7 ! Redistribution and use in source and binary forms, with or without modification, are permitted 8 8 ! provided that the following conditions are met: 9 ! 9 10 10 ! * Redistributions of source code must retain the above copyright notice, this list 11 11 ! of conditions and the following disclaimer. … … 16 16 ! to endorse or promote products derived from this software without specific prior written 17 17 ! permission. 18 ! 18 19 19 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 20 20 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 25 25 ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 26 26 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 27 ! 28 29 ! 27 30 28 ! History: 31 29 ! May 2009 - Robert Pincus - Initial version 32 30 ! Dec 2009 - Robert Pincus - Tiny revisions 33 ! 31 34 32 MODULE MOD_COSP_Modis_Simulator 35 33 USE MOD_COSP_CONSTANTS … … 44 42 !------------------------------------------------------------------------------------------------ 45 43 ! Public type 46 ! 44 47 45 ! Summary statistics from MODIS retrievals 48 46 type COSP_MODIS 49 47 ! Dimensions 50 48 integer :: Npoints ! Number of gridpoints 51 52 ! 49 53 50 ! Grid means; dimension nPoints 54 ! 51 55 52 real, dimension(:), pointer :: & 56 53 Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & … … 61 58 Cloud_Top_Pressure_Total_Mean, & 62 59 Liquid_Water_Path_Mean, Ice_Water_Path_Mean 63 ! 60 64 61 ! Also need the ISCCP-type optical thickness/cloud top pressure histogram 65 ! 62 66 63 real, dimension(:, :, :), pointer :: Optical_Thickness_vs_Cloud_Top_Pressure 67 64 real, dimension(:, :, :), pointer :: Optical_Thickness_vs_ReffICE … … 126 123 integer, dimension(count(gridBox%sunlit(:) <= 0)) :: notSunlit 127 124 ! ------------------------------------------------------------ 128 129 ! 125 130 126 ! Are there any sunlit points? 131 ! 127 132 128 nSunlit = count(gridBox%sunlit(:) > 0) 133 129 if(nSunlit > 0) then … … 135 131 nPoints = gridBox%Npoints 136 132 nSubCols = subCols%Ncolumns 137 ! 133 138 134 ! This is a vector index indicating which points are sunlit 139 ! 135 140 136 sunlit(:) = pack((/ (i, i = 1, nPoints ) /), mask = gridBox%sunlit(:) > 0) 141 137 notSunlit(:) = pack((/ (i, i = 1, nPoints ) /), mask = .not. gridBox%sunlit(:) > 0) 142 143 ! 138 144 139 ! Copy needed quantities, reversing vertical order and removing points with no sunlight 145 ! 140 146 141 pressureLevels(:, 1) = 0.0 ! Top of model, following ISCCP sim 147 142 temperature(:, :) = gridBox%T (sunlit(:), nLevels:1:-1) 148 143 pressureLayers(:, :) = gridBox%p (sunlit(:), nLevels:1:-1) 149 144 pressureLevels(:, 2:) = gridBox%ph(sunlit(:), nLevels:1:-1) 150 151 ! 145 152 146 ! Subcolumn properties - first stratiform cloud... 153 ! 147 154 148 where(subCols%frac_out(sunlit(:), :, :) == I_LSC) 155 149 !opticalThickness(:, :, :) = & … … 180 174 end do 181 175 182 !183 176 ! .. then add convective cloud... 184 ! 177 185 178 where(subCols%frac_out(sunlit(:), :, :) == I_CVC) 186 179 !opticalThickness(:, :, :) = & … … 201 194 end do 202 195 203 !204 196 ! Reverse vertical order 205 ! 197 206 198 opticalThickness(:, :, :) = opticalThickness(:, :, nLevels:1:-1) 207 199 cloudWater (:, :, :) = cloudWater (:, :, nLevels:1:-1) … … 242 234 jointHistogram,jointHistogram2,jointHistogram3) 243 235 ! DJS2015: END 244 245 ! 236 246 237 ! Copy results into COSP structure 247 ! 238 248 239 modisSim%Cloud_Fraction_Total_Mean(sunlit(:)) = cfTotal(:) 249 240 modisSim%Cloud_Fraction_Water_Mean(sunlit(:)) = cfLiquid … … 273 264 modisSim%Optical_Thickness_vs_ReffICE(sunlit(:),2:numModisTauBins+1,:) = jointHistogram2(:, :, :) 274 265 modisSim%Optical_Thickness_vs_ReffLIQ(sunlit(:),2:numModisTauBins+1,:) = jointHistogram3(:, :, :) 275 ! 266 276 267 ! Reorder pressure bins in joint histogram to go from surface to TOA 277 ! 268 278 269 modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:,:,:) = modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, numModisPressureBins:1:-1) 279 270 if(nSunlit < nPoints) then 280 ! 271 281 272 ! Where it's night and we haven't done the retrievals the values are undefined 282 ! 273 283 274 modisSim%Cloud_Fraction_Total_Mean(notSunlit(:)) = R_UNDEF 284 275 modisSim%Cloud_Fraction_Water_Mean(notSunlit(:)) = R_UNDEF … … 310 301 end if 311 302 else 312 ! 303 313 304 ! It's nightime everywhere - everything is undefined 314 ! 305 315 306 modisSim%Cloud_Fraction_Total_Mean(:) = R_UNDEF 316 307 modisSim%Cloud_Fraction_Water_Mean(:) = R_UNDEF … … 350 341 integer, intent(in) :: Npoints ! Number of sampled points 351 342 type(cosp_MODIS), intent(out) :: x 352 ! 343 353 344 ! Allocate minumum storage if simulator not used 354 ! 345 355 346 if (cfg%LMODIS_sim) then 356 347 x%nPoints = nPoints … … 398 389 SUBROUTINE FREE_COSP_MODIS(x) 399 390 type(cosp_MODIS),intent(inout) :: x 400 ! 391 401 392 ! Free space used by cosp_modis variable. 402 ! 403 393 404 394 if(associated(x%Cloud_Fraction_Total_Mean)) deallocate(x%Cloud_Fraction_Total_Mean) 405 395 if(associated(x%Cloud_Fraction_Water_Mean)) deallocate(x%Cloud_Fraction_Water_Mean) … … 439 429 type(cosp_modis), intent(in ) :: orig 440 430 type(cosp_modis), intent( out) :: copy 441 ! 431 442 432 ! Copy a set of grid points from one cosp_modis variable to another. 443 433 ! Should test to be sure ix and iy refer to the same number of grid points 444 ! 434 445 435 integer :: orig_start, orig_end, copy_start, copy_end 446 436 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_radar.F90
r5095 r5099 1 1 ! (c) British Crown Copyright 2008, the Met Office. 2 2 ! All rights reserved. 3 ! 3 4 4 ! Redistribution and use in source and binary forms, with or without modification, are permitted 5 5 ! provided that the following conditions are met: 6 ! 6 7 7 ! * Redistributions of source code must retain the above copyright notice, this list 8 8 ! of conditions and the following disclaimer. … … 13 13 ! to endorse or promote products derived from this software without specific prior written 14 14 ! permission. 15 ! 15 16 16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 17 17 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_simulator.F90
r4619 r5099 3 3 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_simulator.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! to endorse or promote products derived from this software without specific prior written 16 16 ! permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 25 25 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 26 27 !28 27 ! History: 29 28 ! Jul 2007 - A. Bodas-Salcedo - Initial version 30 29 ! Jan 2013 - G. Cesana - Add new variables linked to the lidar cloud phase 31 !32 30 33 31 #include "cosp_defs.h" -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_stats.F90
r4619 r5099 3 3 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_stats.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! to endorse or promote products derived from this software without specific prior written 16 16 ! permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 25 25 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 26 27 !28 27 ! History: 29 28 ! Jul 2007 - A. Bodas-Salcedo - Initial version … … 34 33 ! Jan 2013 - G. Cesana - Added betaperp and temperature arguments 35 34 ! - Added phase 3D/3Dtemperature/Map output variables in diag_lidar 36 ! 37 ! 35 36 38 37 #include "cosp_defs.h" 39 38 MODULE MOD_COSP_STATS -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_types.F90
r5095 r5099 1 1 ! (c) British Crown Copyright 2008, the Met Office. 2 2 ! All rights reserved. 3 ! 3 4 4 ! Redistribution and use in source and binary forms, with or without modification, are permitted 5 5 ! provided that the following conditions are met: 6 ! 6 7 7 ! * Redistributions of source code must retain the above copyright notice, this list 8 8 ! of conditions and the following disclaimer. … … 13 13 ! to endorse or promote products derived from this software without specific prior written 14 14 ! permission. 15 ! 15 16 16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 17 17 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_utils.F90
r3241 r5099 3 3 ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_utils.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! to endorse or promote products derived from this software without specific prior written 16 16 ! permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 25 25 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 26 27 !28 27 ! History: 29 28 ! Jul 2007 - A. Bodas-Salcedo - Initial version 30 !31 29 32 30 MODULE MOD_COSP_UTILS -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_llnl_stats.F90
r5095 r5099 3 3 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/llnl/llnl_stats.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! nor the names of its contributors may be used to endorse or promote products derived from 16 16 ! this software without specific prior written permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 24 24 ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 25 25 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 ! 26 27 27 ! History 28 ! 28 29 29 ! Jan 2013 - G. Cesana - Added betaperp_tot and temp_tot arguments 30 !31 32 30 33 31 MODULE MOD_LLNL_STATS … … 62 60 ! bmin: mimumum value of first bin 63 61 ! bwidth: bin width 64 ! 62 65 63 ! Output: 2D histogram on each horizontal point (Npoints,Nbins,Nlevels) 66 64 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_lmd_ipsl_stats.F90
r5095 r5099 3 3 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/actsim/lmd_ipsl_stats.F90 $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! contributors may be used to endorse or promote products derived from this software without 16 16 ! specific prior written permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 38 38 ,cfad2,srbval,ncat,ntype,lidarcld,lidarcldtype,lidarcldphase,cldlayer & !OPAQ 39 39 ,cldtype,cldlayerphase,lidarcldtmp,parasolrefl,vgrid_z,profSR) !OPAQ !TIBO 40 ! 40 41 41 ! ----------------------------------------------------------------------------------- 42 42 ! Lidar outputs : 43 ! 43 44 44 ! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction) 45 45 ! and phase cloud fraction (3D, low/mid/high/total and 3D temperature) … … 47 47 ! + 48 48 ! Compute CFADs of lidar scattering ratio SR and of depolarization index 49 ! 49 50 50 ! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France). 51 ! 51 52 52 ! December 2008, S. Bony, H. Chepfer and J-L. Dufresne : 53 53 ! - change of the cloud detection threshold S_cld from 3 to 5, for better … … 62 62 ! June 2010, T. Yokohata, T. Nishimura and K. Ogochi 63 63 ! Optimisation of COSP_CFAD_SR 64 ! 64 65 65 ! January 2013, G. Cesana, H. Chepfer: 66 66 ! - Add the perpendicular component of the backscattered signal (pnorm_perp) in the arguments … … 73 73 ! Reference: Cesana G. and H. Chepfer (2013): Evaluation of the cloud water phase 74 74 ! in a climate model using CALIPSO-GOCCP, J. Geophys. Res., doi: 10.1002/jgrd.50376 75 ! 75 76 76 ! ------------------------------------------------------------------------------------ 77 77 … … 133 133 real, parameter :: S_cld_att = 30. ! New threshold for undefine cloud phase detection 134 134 135 136 !137 135 ! c ------------------------------------------------------- 138 136 ! c 0- Initializations 139 137 ! c ------------------------------------------------------- 140 ! 138 141 139 ! Should be modified in future version 142 140 xmax=undef-1.0 … … 174 172 ! c ------------------------------------------------------- 175 173 if (ok_lidar_cfad) then 176 ! 174 177 175 ! c CFADs of subgrid-scale lidar scattering ratios : 178 176 ! c ------------------------------------------------------- … … 224 222 ! S_att: Threshold for full attenuation 225 223 ! S_clr: Threshold for clear-sky layer 226 ! 224 227 225 !--- Input-Outout arguments 228 226 ! x: variable to process (Npoints,Ncolumns,Nlevels), mofified where saturation occurs 229 ! 227 230 228 ! -- Output arguments 231 229 ! srbval : values of the histogram bins … … 528 526 enddo 529 527 endif 530 !531 528 532 529 iz=1 … … 600 597 601 598 !____________________________________________________________________________________________________ 602 ! 599 603 600 !4.1.a Ice: ATBperp above the phase discrimination line 604 601 !____________________________________________________________________________________________________ 605 ! 602 606 603 if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then ! Ice clouds 607 604 ! ICE with temperature above 273,15°K = Liquid (false ice) … … 645 642 646 643 !____________________________________________________________________________________________________ 647 ! 644 648 645 ! 4.1.b Liquid: ATBperp below the phase discrimination line 649 646 !____________________________________________________________________________________________________ 650 ! 647 651 648 else ! Liquid clouds 652 649 ! Liquid with temperature above 231,15°K … … 710 707 ATB(i,ncol,nlev)*epsilon50 + zeta50 711 708 !____________________________________________________________________________________________________ 712 ! 709 713 710 ! 4.2.a Ice: ATBperp above the phase discrimination line 714 711 !____________________________________________________________________________________________________ 715 ! 712 716 713 ! ICE with temperature above 273,15°K = Liquid (false ice) 717 714 if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then ! Ice clouds … … 756 753 757 754 !____________________________________________________________________________________________________ 758 ! 755 759 756 ! 4.2.b Liquid: ATBperp below the phase discrimination line 760 757 !____________________________________________________________________________________________________ 761 ! 758 762 759 else 763 760 ! Liquid with temperature above 231,15°K … … 816 813 817 814 !____________________________________________________________________________________________________ 818 ! 815 819 816 ! Undefined phase: For a cloud located below another cloud with SR>30 820 817 ! see Cesana and Chepfer 2013 Sect.III.2 821 818 !____________________________________________________________________________________________________ 822 ! 819 823 820 if(toplvlsat.ne.0)then 824 821 do nlev=toplvlsat,1,-1 … … 849 846 850 847 !____________________________________________________________________________________________________ 851 ! 848 852 849 ! Computation of final cloud phase diagnosis 853 850 !____________________________________________________________________________________________________ 854 !855 851 856 852 ! Compute the Ice percentage in cloud = ice/(ice+liq) as a function -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_modis_sim.F90
r5095 r5099 4 4 ! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $ 5 5 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/modis_simulator.F90 $ 6 ! 6 7 7 ! Redistribution and use in source and binary forms, with or without modification, are permitted 8 8 ! provided that the following conditions are met: 9 ! 9 10 10 ! * Redistributions of source code must retain the above copyright notice, this list 11 11 ! of conditions and the following disclaimer. … … 16 16 ! to endorse or promote products derived from this software without specific prior written 17 17 ! permission. 18 ! 18 19 19 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 20 20 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND … … 25 25 ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 26 26 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 27 ! 28 29 ! 27 30 28 ! History: 31 29 ! May 2009 - Robert Pincus - Initial version … … 34 32 ! November 2009 - Robert Pincus - Bux fixes and speed-ups after experience with Rick Hemler using AM2 (GFDL) 35 33 ! January 2010 - Robert Pincus - Added high, middle, low cloud fractions 36 ! 37 38 ! 34 39 35 ! Notes on using the MODIS simulator: 40 36 ! *) You may provide either layer-by-layer values of optical thickness at 0.67 and 2.1 microns, or … … 47 43 ! libraries (parameters re_water_min, re_ice_min, etc.) the simulator returns missing values (re_fill) 48 44 49 !50 45 ! When error conditions are encountered this code calls the function complain_and_die, supplied at the 51 46 ! bottom of this module. Users probably want to replace this with something more graceful. 52 ! 47 53 48 module mod_modis_sim 54 49 USE MOD_COSP_TYPES, only: R_UNDEF, ok_debug_cosp … … 56 51 ! ------------------------------ 57 52 ! Algorithmic parameters 58 ! 59 53 60 54 real, parameter :: ice_density = 0.93 ! liquid density is 1. 61 ! 55 62 56 ! Retrieval parameters 63 ! 57 64 58 real, parameter :: min_OpticalThickness = 0.3, & ! Minimum detectable optical thickness 65 59 CO2Slicing_PressureLimit = 700. * 100., & ! Cloud with higher pressures use thermal methods, units Pa … … 73 67 74 68 logical, parameter :: useSimpleReScheme = .false. 75 ! 69 76 70 ! These are the limits of the libraries for the MODIS collection 5 algorithms 77 71 ! They are also the limits used in the fits for g and w0 78 ! 72 79 73 real, parameter :: re_water_min= 4., re_water_max= 30., re_ice_min= 5., re_ice_max= 90. 80 74 integer, parameter :: num_trial_res = 15 ! increase to make the linear pseudo-retrieval of size more accurate … … 82 76 ! logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations? 83 77 ! DJS2015 END 84 ! 78 85 79 ! Precompute near-IR optical params vs size for retrieval scheme 86 ! 80 87 81 integer, private :: i 88 82 real, dimension(num_trial_res), parameter :: & … … 94 88 ! ------------------------------ 95 89 ! Bin boundaries for the joint optical thickness/cloud top pressure histogram 96 ! 90 97 91 integer, parameter :: numTauHistogramBins = 6, numPressureHistogramBins = 7 98 92 … … 103 97 pressureHistogramBoundaries = (/ 0., 18000., 31000., 44000., 56000., 68000., 80000., 1000000. /) 104 98 real, parameter :: highCloudPressureLimit = 440. * 100., lowCloudPressureLimit = 680. * 100. 105 ! 99 106 100 ! For output - nominal bin centers and bin boundaries. On output pressure bins are highest to lowest. 107 ! 101 108 102 integer, private :: k, l 109 103 real, parameter, dimension(2, numTauHistogramBins) :: & … … 155 149 ! subroutine modis_L2_simulator_oneTau, or 156 150 ! 2) Provide ice and liquid optical depths in each layer 157 ! 151 158 152 interface modis_L2_simulator 159 153 module procedure modis_L2_simulator_oneTau, modis_L2_simulator_twoTaus … … 162 156 !------------------------------------------------------------------------------------------------ 163 157 ! MODIS simulator using specified liquid and ice optical thickness in each layer 164 ! 158 165 159 ! Note: this simulator operates on all points; to match MODIS itself night-time 166 160 ! points should be excluded 167 ! 161 168 162 ! Note: the simulator requires as input the optical thickness and cloud top pressure 169 163 ! derived from the ISCCP simulator run with parameter top_height = 1. … … 171 165 ! and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that 172 166 ! alogrithm in this simulator we simply report the values from the ISCCP simulator. 173 ! 167 174 168 subroutine modis_L2_simulator_twoTaus( & 175 169 temp, pressureLayers, pressureLevels, & … … 212 206 nSubcols = size(liquid_opticalThickness, 1) 213 207 nLevels = size(liquid_opticalThickness, 2) 214 215 ! 208 216 209 ! Initial error checks 217 ! 210 218 211 if(any((/ size(ice_opticalThickness, 1), size(waterSize, 1), size(iceSize, 1), & 219 212 size(isccpTau), size(isccpCloudTopPressure), & … … 234 227 235 228 ! --------------------------------------------------- 236 ! 229 237 230 ! Compute the total optical thickness and the proportion due to liquid in each cell 238 ! 231 239 232 where(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) > 0.) 240 233 tauLiquidFraction(:, :) = liquid_opticalThickness(:, :)/(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :)) … … 243 236 end where 244 237 tauTotal(:, :) = liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) 245 246 ! 238 247 239 ! Optical depth retrieval 248 240 ! This is simply a sum over the optical thickness in each layer 249 241 ! It should agree with the ISCCP values after min values have been excluded 250 ! 242 251 243 retrievedTau(:) = sum(tauTotal(:, :), dim = 2) 252 244 253 !254 245 ! Cloud detection - does optical thickness exceed detection threshold? 255 ! 246 256 247 cloudMask = retrievedTau(:) >= min_OpticalThickness 257 258 ! 248 259 249 ! Initialize initial estimates for size retrievals 260 ! 250 261 251 if(any(cloudMask) .and. .not. useSimpleReScheme) then 262 252 g_w(:) = get_g_nir( phaseIsLiquid, trial_re_w(:)) … … 268 258 do i = 1, nSubCols 269 259 if(cloudMask(i)) then 270 ! 260 271 261 ! Cloud top pressure determination 272 262 ! MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal methods for clouds … … 276 266 ! This assumes linear variation in p between levels. Linear in ln(p) is probably better, 277 267 ! though we'd need to deal with the lowest pressure gracefully. 278 ! 268 279 269 retrievedCloudTopPressure(i) = cloud_top_pressure((/ 0., tauTotal(i, :) /), & 280 270 pressureLevels, & 281 271 CO2Slicing_TauLimit) 282 283 284 ! 272 285 273 ! Phase determination - determine fraction of total tau that's liquid 286 274 ! When ice and water contribute about equally to the extinction we can't tell 287 275 ! what the phase is 288 ! 276 289 277 integratedLiquidFraction = weight_by_extinction(tauTotal(i, :), & 290 278 tauLiquidFraction(i, :), & … … 297 285 retrievedPhase(i) = phaseIsUndetermined 298 286 end if 299 300 ! 287 301 288 ! Size determination 302 ! 289 303 290 if(useSimpleReScheme) then 304 291 ! This is the extinction-weighted size considering only the phase we've chosen 305 ! 292 306 293 if(retrievedPhase(i) == phaseIsIce) then 307 294 retrievedSize(i) = weight_by_extinction(ice_opticalThickness(i, :), & … … 323 310 end if 324 311 else 325 ! 312 326 313 ! Values when we don't think there's a cloud. 327 ! 314 328 315 retrievedCloudTopPressure(i) = R_UNDEF 329 316 retrievedPhase(i) = phaseIsNone … … 337 324 ! mimics what MODIS does to first order. 338 325 ! Of course, ISCCP cloud top pressures are in mb. 339 ! 326 340 327 where(cloudMask(:) .and. retrievedCloudTopPressure(:) > CO2Slicing_PressureLimit) & 341 328 retrievedCloudTopPressure(:) = isccpCloudTopPressure * 100. … … 343 330 end subroutine modis_L2_simulator_twoTaus 344 331 !------------------------------------------------------------------------------------------------ 345 ! 332 346 333 ! MODIS simulator: provide a single optical thickness and the cloud ice and liquid contents; 347 334 ! we'll partition this into ice and liquid optical thickness and call the full MODIS simulator 348 ! 335 349 336 subroutine modis_L2_simulator_oneTau( & 350 337 temp, pressureLayers, pressureLevels, & … … 399 386 tauLiquidFraction(:, :) = 0. 400 387 elsewhere 401 ! 388 402 389 ! Geometic optics limit - tau as LWP/re (proportional to LWC/re) 403 ! 390 404 391 ! Modif AI 02 2018 405 392 tauLiquidFraction(:, :) = (cloudWater(:, :)/max(waterSize(:, :), seuil) ) / & … … 682 669 Liquid_Water_Path_Mean, Ice_Water_Path_Mean, & 683 670 Optical_Thickness_vs_Cloud_Top_Pressure) 684 ! 671 685 672 ! Inputs; dimension nPoints, nSubcols 686 ! 673 687 674 integer, dimension(:, :), intent(in) :: phase 688 675 real, dimension(:, :), intent(in) :: cloud_top_pressure, optical_thickness, particle_size 689 ! 676 690 677 ! Outputs; dimension nPoints 691 ! 678 692 679 real, dimension(:), intent(out) :: & 693 680 Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & … … 702 689 ! --------------------------- 703 690 ! Local variables 704 ! 691 705 692 real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units 706 693 integer :: i, j … … 714 701 nPoints = size(phase, 1) 715 702 nSubcols = size(phase, 2) 716 ! 703 717 704 ! Array conformance checks 718 ! 705 719 706 if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1), & 720 707 size(Cloud_Fraction_Total_Mean), size(Cloud_Fraction_Water_Mean), size(Cloud_Fraction_Ice_Mean), & … … 728 715 if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /) /= nSubcols)) & 729 716 call complain_and_die("Some L3 arrays have wrong number of subcolumns") 730 731 732 ! 717 733 718 ! Include only those pixels with successful retrievals in the statistics 734 ! 719 735 720 validRetrievalMask(:, :) = particle_size(:, :) > 0. 736 721 cloudMask = phase(:, :) /= phaseIsNone .and. validRetrievalMask(:, :) 737 722 waterCloudMask = phase(:, :) == phaseIsLiquid .and. validRetrievalMask(:, :) 738 723 iceCloudMask = phase(:, :) == phaseIsIce .and. validRetrievalMask(:, :) 739 ! 724 740 725 ! Use these as pixel counts at first 741 ! 726 742 727 Cloud_Fraction_Total_Mean(:) = real(count(cloudMask, dim = 2)) 743 728 Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2)) … … 747 732 Cloud_Fraction_Low_Mean(:) = real(count(cloudMask .and. cloud_top_pressure > lowCloudPressureLimit, dim = 2)) 748 733 Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:) 749 750 ! 734 751 735 ! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0. 752 ! 736 753 737 where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1. 754 738 where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1. … … 780 764 / Cloud_Fraction_Ice_Mean(:) 781 765 782 !783 766 ! Normalize pixel counts to fraction 784 767 ! The first three cloud fractions have been set to -1 in cloud-free areas, so set those places to 0. 785 ! 768 786 769 Cloud_Fraction_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols) 787 770 Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols) … … 794 777 ! ---- 795 778 ! Joint histogram 796 ! 779 797 780 do i = 1, numTauHistogramBins 798 781 where(cloudMask(:, :)) … … 826 809 real, intent(in) :: tauLimit 827 810 real :: cloud_top_pressure 828 ! 811 829 812 ! Find the extinction-weighted pressure. Assume that pressure varies linearly between 830 813 ! layers and use the trapezoidal rule. 831 ! 832 814 833 815 real :: deltaX, totalTau, totalProduct 834 816 integer :: i … … 839 821 deltaX = tauLimit - totalTau 840 822 totalTau = totalTau + deltaX 841 ! 823 842 824 ! Result for trapezoidal rule when you take less than a full step 843 825 ! tauIncrement is a layer-integrated value 844 ! 826 845 827 totalProduct = totalProduct & 846 828 + pressure(i-1) * deltaX & … … 859 841 real, intent(in) :: tauLimit 860 842 real :: weight_by_extinction 861 ! 843 862 844 ! Find the extinction-weighted value of f(tau), assuming constant f within each layer 863 ! 864 845 865 846 real :: deltaX, totalTau, totalProduct 866 847 integer :: i … … 892 873 ice_g(:) = get_g_nir( phaseIsIce, ice_size) 893 874 ice_w0(:) = get_ssa_nir(phaseIsIce, ice_size) 894 ! 875 895 876 ! Combine ice and water optical properties 896 ! 877 897 878 g(:) = 0; w0(:) = 0. 898 879 tau(:) = ice_tau(:) + water_tau(:) … … 913 894 real, intent(in) :: tau, obs_Refl_nir 914 895 real :: retrieve_re 915 ! 896 916 897 ! Finds the re that produces the minimum mis-match between predicted and observed reflectance in 917 898 ! MODIS band 7 (near IR) … … 920 901 ! two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0 921 902 ! adding-doubling for total reflectance 922 ! 923 ! 924 ! 903 904 905 925 906 ! Local variables 926 ! 907 927 908 real, parameter :: min_distance_to_boundary = 0.01 928 909 real :: re_min, re_max, delta_re … … 946 927 w0(:) = w0_i(:) 947 928 end if 948 ! 929 949 930 ! 1st attempt at index: w/coarse re resolution 950 ! 931 951 932 predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:)) 952 933 retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) 953 ! 934 954 935 ! If first retrieval works, can try 2nd iteration using greater re resolution 955 ! 936 956 937 ! DJS2015: Remove unused piece of code 957 938 ! if(use_two_re_iterations .and. retrieve_re > 0.) then … … 959 940 ! re_max = retrieve_re + delta_re 960 941 ! delta_re = (re_max - re_min)/real(num_trial_res-1) 961 ! 942 962 943 ! trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /) 963 944 ! g(:) = get_g_nir( phase, trial_re(:)) … … 977 958 real, intent(in) :: yobs 978 959 real :: interpolate_to_min 979 ! 960 980 961 ! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs) 981 962 ! y must be monotonic in x 982 ! 963 983 964 real, dimension(size(x)) :: diff 984 965 integer :: nPoints, minDiffLoc, lowerBound, upperBound … … 1005 986 1006 987 if(diff(lowerBound) * diff(upperBound) < 0) then 1007 ! 988 1008 989 ! Interpolate the root position linearly if we bracket the root 1009 ! 990 1010 991 interpolate_to_min = x(upperBound) - & 1011 992 diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound)) … … 1020 1001 ! -------------------------------------------- 1021 1002 elemental function get_g_nir (phase, re) 1022 ! 1003 1023 1004 ! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function 1024 1005 ! of size for ice and water 1025 1006 ! Fits from Steve Platnick 1026 !1027 1007 1028 1008 integer, intent(in) :: phase … … 1056 1036 real, intent(in) :: re 1057 1037 real :: get_ssa_nir 1058 ! 1038 1059 1039 ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function 1060 1040 ! of size for ice and water 1061 1041 ! Fits from Steve Platnick 1062 ! 1042 1063 1043 real, dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/) 1064 1044 real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /) … … 1107 1087 integer :: i 1108 1088 ! --------------------------------------- 1109 ! 1089 1110 1090 ! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation 1111 ! 1091 1112 1092 cloudMask = tau(:) > 0. 1113 1093 cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask) … … 1125 1105 real, intent(in) :: tauint, gint, w0int 1126 1106 real, intent(out) :: ref, tra 1127 ! 1107 1128 1108 ! Compute reflectance in a single layer using the two stream approximation 1129 1109 ! The code itself is from Lazaros Oreopoulos via Steve Platnick 1130 ! 1110 1131 1111 ! ------------------------ 1132 1112 ! Local variables … … 1137 1117 real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & 1138 1118 rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th 1139 ! 1119 1140 1120 ! Compute reflectance and transmittance in a single layer using the two stream approximation 1141 1121 ! The code itself is from Lazaros Oreopoulos via Steve Platnick 1142 ! 1122 1143 1123 f = gint**2 1144 1124 tau = (1 - w0int * f) * tauint … … 1204 1184 real, intent(in) :: tauint, gint, w0int 1205 1185 real :: two_stream_reflectance 1206 ! 1186 1207 1187 ! Compute reflectance in a single layer using the two stream approximation 1208 1188 ! The code itself is from Lazaros Oreopoulos via Steve Platnick 1209 ! 1189 1210 1190 ! ------------------------ 1211 1191 ! Local variables … … 1279 1259 real, dimension(:), intent(in) :: Refl, Tran 1280 1260 real, intent(out) :: Refl_tot, Tran_tot 1281 ! 1261 1282 1262 ! Use adding/doubling formulas to compute total reflectance and transmittance from layer values 1283 ! 1284 1263 1285 1264 integer :: i 1286 1265 real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/optics_lib.F90
r5095 r5099 4 4 ! Compiled/Modified: 5 5 ! 07/01/06 John Haynes (haynes@atmos.colostate.edu) 6 ! 6 7 7 ! m_wat (subroutine) 8 8 ! m_ice (subroutine) … … 19 19 subroutine m_wat(freq, tk, n_r, n_i) 20 20 implicit none 21 ! 21 22 22 ! Purpose: 23 23 ! compute complex index of refraction of liquid water 24 ! 24 25 25 ! Inputs: 26 26 ! [freq] frequency (GHz) 27 27 ! [tk] temperature (K) 28 ! 28 29 29 ! Outputs: 30 30 ! [n_r] real part index of refraction 31 31 ! [n_i] imaginary part index of refraction 32 ! 32 33 33 ! Reference: 34 34 ! Based on the work of Ray (1972) 35 ! 35 36 36 ! Coded: 37 37 ! 03/22/05 John Haynes (haynes@atmos.colostate.edu) … … 83 83 subroutine m_ice(freq,t,n_r,n_i) 84 84 implicit none 85 ! 85 86 86 ! Purpose: 87 87 ! compute complex index of refraction of ice 88 ! 88 89 89 ! Inputs: 90 90 ! [freq] frequency (GHz) 91 91 ! [t] temperature (K) 92 ! 92 93 93 ! Outputs: 94 94 ! [n_r] real part index of refraction 95 95 ! [n_i] imaginary part index of refraction 96 ! 96 97 97 ! Reference: 98 98 ! Fortran 90 port from IDL of REFICE by Stephen G. Warren 99 ! 99 100 100 ! Modified: 101 101 ! 05/31/05 John Haynes (haynes@atmos.colostate.edu) … … 121 121 ! Allowable wavelength range extends from 0.045 microns to 8.6 meter 122 122 ! temperature dependence only considered beyond 167 microns. 123 ! 123 124 124 ! interpolation is done n_r vs. log(xlam) 125 125 ! n_r vs. t 126 126 ! log(n_i) vs. log(xlam) 127 127 ! log(n_i) vs. t 128 ! 128 129 129 ! Stephen G. Warren - 1983 130 130 ! Dept. of Atmospheric Sciences 131 131 ! University of Washington 132 132 ! Seattle, Wa 98195 133 ! 133 134 134 ! Based on 135 ! 135 136 136 ! Warren,S.G.,1984. 137 137 ! Optical constants of ice from the ultraviolet to the microwave. 138 138 ! Applied Optics,23,1206-1225 139 ! 139 140 140 ! Reference temperatures are -1.0,-5.0,-20.0, and -60.0 deg C 141 141 … … 576 576 ! subroutine MIEINT 577 577 ! ---------------------------------------------------------------------------- 578 ! 578 579 579 ! General purpose Mie scattering routine for single particles 580 580 ! Author: R Grainger 1990 … … 583 583 ! code to ensure correct calculation of backscatter coeficient 584 584 ! Options/Extend_Source 585 ! 585 586 586 Subroutine MieInt(Dx, SCm, Inp, Dqv, Dqxt, Dqsc, Dbsc, Dg, Xs1, Xs2, DPh, Error) 587 587 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/pf_to_mr.F
r5095 r5099 3 3 ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/llnl/pf_to_mr.f $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! nor the names of its contributors may be used to endorse or promote products derived from 16 16 ! this software without specific prior written permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/prec_scops.F
r5095 r5099 3 3 ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ 4 4 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/llnl/prec_scops.f $ 5 ! 5 6 6 ! Redistribution and use in source and binary forms, with or without modification, are permitted 7 7 ! provided that the following conditions are met: 8 ! 8 9 9 ! * Redistributions of source code must retain the above copyright notice, this list 10 10 ! of conditions and the following disclaimer. … … 15 15 ! nor the names of its contributors may be used to endorse or promote products derived from 16 16 ! this software without specific prior written permission. 17 ! 17 18 18 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 19 19 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator.F90
r5095 r5099 17 17 18 18 ! Purpose: 19 ! 19 20 20 ! Simulates a vertical profile of radar reflectivity 21 21 ! Originally Part of QuickBeam v1.04 by John Haynes & Roger Marchand. 22 22 ! but has been substantially modified since that time by 23 23 ! Laura Fowler and Roger Marchand (see modifications below). 24 ! 24 25 25 ! Inputs: 26 ! 26 27 27 ! [hp] structure that defines hydrometeor types and other radar properties 28 ! 28 29 29 ! [nprof] number of hydrometeor profiles 30 30 ! [ngate] number of vertical layers 31 ! 31 32 32 ! [undef] missing data value 33 33 ! (The following 5 arrays must be in order from closest to the radar 34 34 ! to farthest...) 35 ! 35 36 36 ! [hgt_matrix] height of hydrometeors (km) 37 37 ! [p_matrix] pressure profile (hPa) 38 38 ! [t_matrix] temperature profile (K) 39 39 ! [rh_matrix] relative humidity profile (%) -- only needed if gaseous aborption calculated. 40 ! 40 41 41 ! [hm_matrix] table of hydrometeor mixing rations (g/kg) 42 42 ! [re_matrix] table of hydrometeor effective radii. 0 ==> use defaults. (units=microns) 43 43 ! [Np_matrix] table of hydrometeor number concentration. 0 ==> use defaults. (units = 1/kg) 44 ! 44 45 45 ! Outputs: 46 ! 46 47 47 ! [Ze_non] radar reflectivity without attenuation (dBZ) 48 48 ! [Ze_ray] Rayleigh reflectivity (dBZ) … … 50 50 ! [g_atten_to_vol] gaseous atteunation, radar to vol (dB) 51 51 ! [dBZe] effective radar reflectivity factor (dBZ) 52 ! 52 53 53 ! Optional: 54 54 ! [g_to_vol_in] integrated atten due to gases, r>v (dB). … … 57 57 ! [g_to_vol_out] integrated atten due to gases, r>v (dB). 58 58 ! If present then gaseous absorption for each profile is returned here. 59 ! 59 60 60 ! Created: 61 61 ! 11/28/2005 John Haynes (haynes@atmos.colostate.edu) 62 ! 62 63 63 ! Modified: 64 64 ! 09/2006 placed into subroutine form (Roger Marchand,JMH) … … 66 66 ! 01/2008 'Do while' to determine if hydrometeor(s) present in volume 67 67 ! changed for vectorization purposes (A. Bodas-Salcedo) 68 ! 68 69 69 ! 07/2010 V3.0 ... Modified to load or save scale factors to disk as a Look-Up Table (LUT) 70 70 ! ... All hydrometeor and radar simulator properties now included in hp structure … … 73 73 ! Also ... Support of Morrison 2-moment style microphyscis (Np_matrix) added 74 74 ! ... Changes implement by Roj Marchand following work by Laura Fowler 75 ! 75 76 76 ! 10/2011 Modified ngate loop to go in either direction depending on flag 77 77 ! hp%radar_at_layer_one. This affects the direction in which attenuation is summed. 78 ! 78 79 79 ! Also removed called to AVINT for gas and hydrometeor attenuation and replaced with simple 80 80 ! summation. (Roger Marchand) 81 ! 82 ! 81 82 83 83 ! ----- INPUTS ----- 84 84 … … 153 153 g_to_vol_out_present = present(g_to_vol_out) 154 154 155 !156 155 ! load scaling matricies from disk -- but only the first time this subroutine is called 157 ! 156 158 157 if(hp%load_scale_LUTs) then 159 158 call load_scale_LUTs(hp) … … 223 222 if(Re.gt.0) then 224 223 ! determine index in to scale LUT 225 ! 224 226 225 ! distance between Re points (defined by "base" and "step") for 227 226 ! each interval of size Re_BIN_LENGTH … … 278 277 ! alternative is to comment out above two lines and use the following block 279 278 ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density 280 ! 279 281 280 ! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3) !MG Mie approach 282 281 … … 373 372 374 373 ! :: attenuation due to hydrometeors between radar and volume 375 ! 374 376 375 ! NOTE old scheme integrates attenuation only for the layers ABOVE 377 376 ! the current layer ... i.e. 1 to k-1 rather than 1 to k ... -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator_init.F90
r5095 r5099 12 12 13 13 ! Purpose: 14 ! 14 15 15 ! Initialize variables used by the radar simulator. 16 16 ! Part of QuickBeam v3.0 by John Haynes and Roj Marchand 17 ! 18 ! 17 18 19 19 ! Inputs: 20 20 ! [] from data in hydrometeor class input 21 ! 21 22 22 ! [freq] radar frequency (GHz) 23 ! 23 24 24 ! [k2] |K|^2, the dielectric constant, set to -1 to use the 25 25 ! frequency dependent default 26 ! 26 27 27 ! [use_gas_abs] 1=do gaseous abs calcs, 0=no gasesous absorbtion calculated, 28 28 ! 2=calculate (and use) absorption for first profile on all profiles 29 ! 29 30 30 ! [undef] mising data value 31 31 ! [nhclass] number of hydrometeor types 32 ! 32 33 33 ! For each hydrometero type: 34 34 ! hclass_type Type of distribution (see quickbeam documentation) 35 35 ! hclass_phase 1==ice, 0=liquid 36 ! 36 37 37 ! hclass_dmin minimum diameter allowed is drop size distribution N(D<Dmin)=0 38 38 ! hclass_dmax maximum diameter allowed is drop size distribution N(D>Dmax)=0 39 ! 39 40 40 ! hclass_apm,hclass_bpm Density of partical apm*D^bpm or constant = rho 41 41 ! hclass_rho, 42 ! 42 43 43 ! hclass_p1,hclass_p2, Default values of DSD parameters (see quickbeam documentation) 44 44 ! hclass_p3, 45 ! 45 46 46 ! load_scale_LUTs_flag Flag, load scale factors Look Up Table from file at start of run 47 47 ! update_scale_LUTs_flag Flag, save new scale factors calculated during this run to LUT 48 48 ! LUT_file_name Name of file containing LUT 49 ! 49 50 50 ! Outputs: 51 51 ! [hp] structure that define hydrometeor types 52 ! 52 53 53 ! Modified: 54 54 ! 08/23/2006 placed into subroutine form (Roger Marchand) … … 75 75 integer :: i,j 76 76 real*8 :: delt, deltp 77 78 ! 77 79 78 ! set radar simulation properites 80 ! 79 81 80 hp%freq=freq 82 81 hp%k2=k2 … … 88 87 hp%update_scale_LUTs=update_scale_LUTs_flag 89 88 hp%scale_LUT_file_name=LUT_file_name 90 91 ! 89 92 90 ! Store settings for hydrometeor types in hp (class_parameter) structure. 93 ! 91 94 92 do i = 1,nhclass 95 93 hp%dtype(i) = hclass_type(i) … … 104 102 hp%p3(i) = hclass_p3(i) 105 103 enddo 106 107 ! 104 108 105 ! initialize scaling array 109 ! 106 110 107 hp%N_scale_flag = .false. 111 108 hp%fc = undef … … 116 113 hp%Zr_scaled = 0.0 117 114 hp%kr_scaled = 0.0 118 119 ! 115 120 116 ! set up Re bin "structure" for z_scaling 121 ! 117 122 118 hp%base_list(1)=0; 123 119 do j=1,Re_MAX_BIN … … 131 127 enddo 132 128 133 !134 129 ! set up Temperature bin structure used for z scaling 135 ! 130 136 131 do i=1,cnt_ice 137 132 hp%mt_tti(i)=(i-1)*5-90 + 273.15 … … 141 136 hp%mt_ttl(i)=(i-1)*5-60 + 273.15 142 137 enddo 143 144 ! 138 145 139 ! define array discrete diameters used in mie calculations 146 ! 140 147 141 delt = (log(dmax)-log(dmin))/(nd-1) 148 142 deltp = exp(delt) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/scale_luts_io.F90
r3233 r5099 1 1 ! scale_LUT_io: Contains subroutines to load and save scaling Look Up Tables (LUTs) to a file 2 ! 2 3 3 ! June 2010 Written by Roj Marchand 4 4 … … 16 16 logical :: LUT_file_exists 17 17 integer :: i,j,k,ind 18 19 ! 18 20 19 ! load scale LUT from file 21 ! 20 22 21 inquire(file=trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', & 23 22 exist=LUT_file_exists) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/scops.F
r5095 r5099 8 8 ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ 9 9 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/icarus-scops-4.1-bsd/scops.f $ 10 ! 10 11 11 ! Redistribution and use in source and binary forms, with or without 12 12 ! modification, are permitted provided that the 13 13 ! following conditions are met: 14 ! 14 15 15 ! * Redistributions of source code must retain the above 16 16 ! copyright notice, this list of conditions and the following … … 24 24 ! derived from this software without specific prior written 25 25 ! permission. 26 ! 26 27 27 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 28 28 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT … … 36 36 ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 37 37 ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 38 ! 38 39 39 ! *****************************COPYRIGHT******************************* 40 40 ! *****************************COPYRIGHT******************************* -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/zeff.F90
r5095 r5099 10 10 ! Part of QuickBeam v1.03 by John Haynes 11 11 ! http://reef.atmos.colostate.edu/haynes/radarsim 12 ! 12 13 13 ! Inputs: 14 14 ! [freq] radar frequency (GHz) … … 23 23 ! [qs] ... efficiencies; otherwise set to -1 24 24 ! [rho_e] medium effective density (kg m^-3) (-1 = pure) 25 ! 25 26 26 ! Outputs: 27 27 ! [z_eff] unattenuated effective reflectivity factor (mm^6/m^3) 28 28 ! [z_ray] reflectivity factor, Rayleigh only (mm^6/m^3) 29 29 ! [kr] attenuation coefficient (db km^-1) 30 ! 30 31 31 ! Created: 32 32 ! 11/28/05 John Haynes (haynes@atmos.colostate.edu)
Note: See TracChangeset
for help on using the changeset viewer.