source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 1292

Last change on this file since 1292 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

  • Optimized computations in paramfoto_compact (twice less dlog10 calculations)
  • Checked consistency before/after modification in debug mode
  • Checked performance is not impacted (same as before)
File size: 20.0 KB
Line 
1      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
2   ! to use  'getin'
3      USE ioipsl_getincom
4      use netcdf
5      use tracer_mod, only: nqmx, noms
6      use comgeomfi_h, only: long, lati
7      use surfdat_h, only: watercaptag, frost_albedo_threshold,
8     &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
9     &                     albedice, dryness
10#ifndef MESOSCALE
11      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
12      use mod_phys_lmdz_para, only : is_master, gather, scatter
13#endif
14      USE comcstfi_h
15      IMPLICIT NONE
16c=======================================================================
17c
18c   creation des calottes pour l'etat initial
19c
20c=======================================================================
21c-----------------------------------------------------------------------
22c   Declarations:
23c   -------------
24#include "dimensions.h"
25#include "callkeys.h"
26#include "datafile.h"
27
28      integer,intent(in) :: ngrid ! number of atmospheric columns
29      real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
30      real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
31      real,intent(out) :: psolaralb(ngrid,2) ! albedo
32
33      INTEGER ig,icap,iq,alternate
34      REAL icedryness ! ice dryness
35     
36      ! longwatercaptag is watercaptag. Trick for some compilers
37      LOGICAL, DIMENSION(100000) :: longwatercaptag
38     
39! There are 3 different modes for ice distribution:
40! icelocationmode = 1 ---> based on data from surface.nc
41! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
42! icelocationmode = 3 ---> based on logical relations for latitude and longitude
43! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
44      INTEGER,SAVE :: icelocationmode = 2
45       
46       
47      !in case icelocationmode == 1
48      INTEGER i,j
49      INTEGER     imd,jmd
50      PARAMETER   (imd=360,jmd=180)
51      REAL        zdata(imd,jmd)
52      REAL        zelat,zelon
53
54#ifndef MESOSCALE
55      INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
56#endif
57      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
58
59      REAL step,count,ratiolat
60
61      INTEGER   ierr,nid,nvarid
62     
63      REAL,SAVE :: min_icevalue = 500.
64      character(len=50) :: string = 'thermal'
65     
66      character (len=100) :: zedatafile
67
68#ifdef MESOSCALE
69
70      do ig=1,ngrid
71
72         !write(*,*) "all qsurf to zero. dirty."
73         do iq=1,nqmx
74         qsurf(ig,iq)=0.  !! on jette les inputs GCM
75                          !! on regle juste watercaptag
76                          !! il faudrait garder les inputs GCM
77                          !! si elles sont consequentes
78         enddo
79         if ( ( lati(ig)*180./pi .gt. 70. ) .and.
80     .        ( albedodat(ig) .ge. 0.26   ) )  then
81                 write(*,*)"outlier ",ig
82                 watercaptag(ig)  = .true.
83                 dryness(ig)      = 1.
84                 albedodat(ig)    = albedo_h2o_ice  !! pour output
85         else
86                 watercaptag(ig)  = .false.
87                 dryness(ig)      = 1.
88         endif
89
90      enddo
91#endif
92! problem with nested precompiling flags
93
94#ifndef MESOSCALE
95      ! to handle parallel cases
96#if CPP_PARA
97      logical watercaptag_glo(klon_glo)
98      real dryness_glo(klon_glo)
99      real lati_glo(klon_glo)
100      real long_glo(klon_glo)
101#else
102      logical watercaptag_glo(ngrid)
103      real dryness_glo(ngrid)
104      real lati_glo(ngrid)
105      real long_glo(ngrid)
106#endif
107#endif
108
109#ifndef MESOSCALE
110
111c
112c=======================================================================
113! Initialize watercaptag (default is false)
114      watercaptag_glo(:)=.false.
115
116c     water ice outliers
117c     ------------------------------------------
118
119      IF ((water) .and. (caps)) THEN
120     
121c Perennial H20 north cap defined by watercaptag=true (allows surface to be
122c hollowed by sublimation in vdifc).
123
124c We might not want albedodat to be modified because it is used to write
125c restart files. Instead, albedo is directly modified when needed (i.e.
126c if we have watercaptag and no co2 ice), below and in albedocaps.F90
127
128c       "Dryness coefficient" controlling the evaporation and
129c        sublimation from the ground water ice (close to 1)
130c        HERE, the goal is to correct for the fact
131c        that the simulated permanent water ice polar caps
132c        is larger than the actual cap and the atmospheric
133c        opacity not always realistic.
134
135         alternate = 0
136         
137         if (ngrid .ne. 1) then
138           watercaptag(:) = .false.
139           longwatercaptag(:) = .false.
140         endif
141         
142         write(*,*) "surfini: Ice dryness ?"
143         icedryness=1. ! default value
144         call getin("icedryness",icedryness)
145         write(*,*) "surfini: icedryness = ",icedryness
146         dryness (:) = icedryness
147         
148      ! To be able to run in parallel, we work on the full grid
149      ! and dispatch results afterwards
150
151      ! start by geting latitudes and logitudes on full grid
152      ! (in serial mode, this is just a copy)
153      call gather(lati,lati_glo)
154      call gather(long,long_glo)
155
156      if (is_master) then
157
158        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
159     
160         print*, 'ngrid = 1, do no put ice caps in surfini.F'
161
162        ELSE IF (icelocationmode .eq. 1) THEN
163     
164         print*,'Surfini: ice caps defined from surface.nc'
165           
166! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
167! Typically, it is for thermal inertia above 500 tiu.
168! Two conditions are verified:
169! 1. GCM ice caps are defined such as area is conserved for a given latitude
170! (the approximation is that all points within the GCM latitude resolution have the same area).
171! 2. caps are placed to fill the GCM points with the most detected ice first.
172     
173
174           
175         zedatafile = trim(datafile)
176 
177       
178         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
179     &   NF90_NOWRITE,nid)
180     
181         IF (ierr.NE.nf90_noerr) THEN
182       write(*,*)'Error : cannot open file surface.nc '
183       write(*,*)'(in phymars/surfini.F)'
184       write(*,*)'It should be in :',trim(zedatafile),'/'
185       write(*,*)'1) You can set this path in the callphys.def file:'
186       write(*,*)'   datadir=/path/to/the/datafiles'
187       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
188       write(*,*)'   can be obtained online on:'
189       write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
190       CALL ABORT
191         ENDIF
192     
193     
194         ierr=nf90_inq_varid(nid, string, nvarid)
195         if (ierr.ne.nf90_noerr) then
196          write(*,*) 'surfini error, cannot find ',trim(string)
197          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
198          write(*,*)trim(nf90_strerror(ierr))
199          stop
200         endif
201
202         ierr=nf90_get_var(nid, nvarid, zdata)
203
204         if (ierr.ne.nf90_noerr) then
205          write(*,*) 'surfini: error failed loading ',trim(string)
206          write(*,*)trim(nf90_strerror(ierr))
207          stop
208         endif
209 
210                     
211         ierr=nf90_close(nid)
212 
213
214         nb_ice(:,1) = 1 ! default: there is no ice
215         latice(:,1) = 1
216         lonice(:,1) = 1
217         nb_ice(:,2) = 0
218         latice(:,2) = 0
219         lonice(:,2) = 0
220         !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
221
222         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
223         do ig=2,klon_glo-1
224     
225          ! loop over the surface file grid     
226          do i=1,imd
227           do j=1,jmd
228             zelon = i - 180.
229             zelat = 90. - j
230            if ((abs(lati_glo(ig)*180./pi-zelat).le.90./real(jjm)) .and.
231     &        (abs(long_glo(ig)*180./pi-zelon).le.180./real(iim))) then
232              ! count all points in that GCM grid point
233              nb_ice(ig,1) = nb_ice(ig,1) + 1
234              if (zdata(i,j) > min_icevalue)
235                 ! count all detected points in that GCM grid point
236     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
237             endif
238           enddo
239          enddo 
240
241        ! projection of nb_ice on GCM lat and lon axes
242          latice(1+(ig-2)/iim,:) =
243     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
244          lonice(1+mod(ig-2,iim),:) =
245     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
246
247         enddo ! of do ig=2,klon_glo-1
248     
249
250     
251         ! special case for poles
252         nb_ice(1,2)   = 1  ! ice prescribed on north pole
253         latice(1,:)   = nb_ice(1,:)
254         lonice(1,:)   = nb_ice(1,:)
255         latice(jjm,:) = nb_ice(ngrid,:)
256         lonice(iim,:) = nb_ice(ngrid,:)
257     
258     
259!      print*, 'latice TOT', latice(:,1)
260!      print*, 'latice FOUND', latice(:,2)
261!      print*, 'lonice TOT', lonice(:,1)
262!      print*, 'lonice FOUND', lonice(:,2)
263     
264!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
265!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
266     
267!      print*,''
268!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
269!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
270     
271   
272         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
273         do i=1,jjm/2
274          step  = 1. ! threshold to add ice cap
275          count = 0. ! number of ice GCM caps at this latitude
276          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
277          ratiolat  = real(latice(i,2))/real(latice(i,1))
278          !print*,'i',i,(i-1)*iim+2,i*iim+1
279     
280          ! put ice caps while there is not enough ice,
281          ! as long as the threshold is above 20%
282          do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
283           count = 0.
284           ! loop over GCM longitudes
285           do j=1,iim
286            ! if the detected ice ratio in the GCM grid point
287            ! is more than 'step', then add ice
288            if (real(nb_ice((i-1)*iim+1+j,2))
289     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
290                  watercaptag_glo((i-1)*iim+1+j) = .true.
291                  count = count + 1
292            endif
293           enddo ! of do j=1,iim
294           !print*, 'step',step,count,ratiolat*iim
295           step = step - 0.01
296          enddo ! of do while
297          !print*, 'step',step,count,ratiolat*iim
298
299         enddo ! of do i=1,jjm/2
300           
301
302        ELSE IF (icelocationmode .eq. 2) THEN
303     
304         print*,'Surfini: predefined ice caps'
305     
306         if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
307           
308          print*,'water ice caps distribution for 32x24 resolution'
309          longwatercaptag(1:9)    = .true. ! central cap - core
310          longwatercaptag(26:33)  = .true. ! central cap
311          longwatercaptag(1:33)  = .true. ! central cap
312          longwatercaptag(56)  = .true. ! central cap
313          longwatercaptag(58)  = .true. ! central cap
314          longwatercaptag(60)  = .true. ! central cap
315          longwatercaptag(62)  = .true. ! central cap
316          longwatercaptag(64)  = .true. ! central cap
317!---------------------   OUTLIERS  ----------------------------
318
319         else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
320
321          print*,'water ice caps distribution for 64x48 resolution'
322          longwatercaptag(1:65)   = .true. ! central cap - core
323          longwatercaptag(75:85)  = .true. ! central cap
324          longwatercaptag(93:114) = .true. ! central cap
325!---------------------   OUTLIERS  ----------------------------
326          if (.true.) then
327          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
328          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
329          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
330          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
331          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
332          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
333          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
334          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
335          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
336          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
337          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
338          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
339          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
340          longwatercaptag(194)    = .true. ! outlier, lat = 75
341          longwatercaptag(203)    = .true. ! outlier, lat = 75
342          longwatercaptag(207)    = .true. ! outlier, lat = 75
343          longwatercaptag(244)    = .true. ! outlier, lat = 75
344          longwatercaptag(246)    = .true. ! outlier, lat = 75
345          longwatercaptag(250)    = .true. ! outlier, lat = 75
346          longwatercaptag(252)    = .true. ! outlier, lat = 75
347          longwatercaptag(254)    = .true. ! outlier, lat = 75
348          longwatercaptag(256)    = .true. ! outlier, lat = 75
349          endif
350!--------------------------------------------------------------       
351
352
353           
354         else if (klon_glo .ne. 1) then
355       
356          print*,'No predefined ice location for this resolution :',
357     &           iim,jjm
358          print*,'Please change icelocationmode in surfini.F'
359          print*,'Or add some new definitions ...'
360          call abort
361         
362         endif
363
364         do ig=1,klon_glo
365          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
366         enddo
367
368
369        ELSE IF (icelocationmode .eq. 3) THEN
370     
371         print*,'Surfini: ice caps defined by lat and lon values'
372
373         do ig=1,klon_glo
374         
375c-------- Towards olympia planitia water caps -----------
376c--------------------------------------------------------
377
378          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
379     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
380     .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
381     .           ( long_glo(ig)*180./pi .le. 181. ) )
382     .         .or.
383
384     .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
385     .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
386     .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
387     .           ( long_glo(ig)*180./pi .le. 168. ) )
388     .         .or.
389     .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
390     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
391     .           ( long_glo(ig)*180./pi .ge. -150.) .and.
392     .           ( long_glo(ig)*180./pi .le. -110.) ) )
393     .         then
394             
395               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
396              !    watercaptag(ig)=.true.
397                  alternate = 1
398               else
399                  alternate = 0
400               endif !end if alternate = 0
401               
402          endif
403
404c----------- Opposite olympia planitia water cap --------
405c--------------------------------------------------------
406
407          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
408     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
409     .         .and.
410     .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
411     .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
412!     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
413!     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
414!     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
415!     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
416        !   watercaptag_glo(ig)=.true.
417          endif
418
419
420c -------------------- Central cap ----------------------
421c--------------------------------------------------------
422
423          if (abs(lati_glo(ig)*180./pi).gt.80)
424     .          watercaptag_glo(ig)=.true.
425           
426c--------------------------------------------------------
427c--------------------------------------------------------
428         end do ! of (klon_glo)
429
430
431        ELSE
432     
433         print*, 'In surfini.F, icelocationmode is ', icelocationmode
434         print*, 'It should be 1, 2 or 3.'
435         call abort
436
437        ENDIF ! of if (icelocation)
438       
439       
440        ! print caps locations - useful for plots too
441        print*,'surfini: latitude | longitude | ig'
442        do ig=1,klon_glo
443          dryness_glo(ig) = icedryness
444
445          if (watercaptag_glo(ig)) then
446             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
447     &              long_glo(ig)*180./pi, ig
448          endif
449        enddo
450       
451       endif !of if (is_master)
452       
453       ! Now scatter fields watercaptag and dryness from master to all
454       ! (is just a plain copy in serial mode)
455       call scatter(dryness_glo,dryness)
456       call scatter(watercaptag_glo,watercaptag)
457       
458! end of #else of #ifdef MESOSCALE
459       ENDIF ! (caps & water)
460#endif       
461
462c ===============================================================
463c      INITIAL ALBEDO
464c ===============================================================
465
466         write(*,*)"surfini: water frost thickness",
467     s     frost_albedo_threshold
468         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
469         write(*,*)"surfini: water ice TI:", inert_h2o_ice
470
471c        To start with : Initial albedo = observed dataset
472c        -------------------------------------------------
473         DO ig=1,ngrid
474              psolaralb(ig,1)=albedodat(ig)
475              psolaralb(ig,2)=albedodat(ig)
476         END DO
477         PRINT*,'surfini: minimum albedo without water caps',
478     &          minval(albedodat)
479         PRINT*,'surfini: maximum albedo without water caps',
480     &          maxval(albedodat)
481
482c        initial albedo if permanent H2O ice is present
483c        ------------------------------------------------
484         IF ((water) .and. (caps)) THEN
485           DO ig=1,ngrid
486            IF (watercaptag(ig)) THEN
487              psolaralb(ig,1) = albedo_h2o_ice
488              psolaralb(ig,2) = albedo_h2o_ice
489            ENDIF
490           END DO
491           PRINT*,'surfini: minimum albedo with water caps',
492     &            minval(albedodat)
493           PRINT*,'surfini: maximum albedo with water caps',
494     &            maxval(albedodat)
495
496         ENDIF
497
498c      changing initial albedo if CO2 ice is present
499c      -------------------------------------------
500
501       DO ig=1,ngrid
502         IF (piceco2(ig) .GT. 0.) THEN
503             IF(lati(ig).LT. 0.) THEN
504                icap=2 ! Southern hemisphere
505             ELSE
506                icap=1 ! Northern hemisphere
507             ENDIF
508             psolaralb(ig,1) = albedice(icap)
509             psolaralb(ig,2) = albedice(icap)   
510         END IF
511       END DO
512
513c      changing initial albedo if water ice frost is present
514c      -------------------------------------------
515       IF (water) THEN
516          do iq=1,nqmx
517c          if there is frost and surface albedo is set to albedo_h2o_ice
518           if(noms(iq).eq."h2o_ice") then
519             do ig=1,ngrid
520             
521              if ((watercaptag(ig).eqv..false.)
522     &     .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
523              print*, ''
524              print*, '!!! PROBLEM in SURFINI !!!!'
525              print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
526     & WATERCAPTAG IS FALSE'
527              print*, ''
528              print*, 'ig,qsurf,threshold' ,
529     &         ig, qsurf(ig,iq), -frost_albedo_threshold
530              print*, ''
531              print*, '1) Check h2o_ice in startfi and ice
532     & distribution in surfini'
533              print*, '2) Use ini_h2osurf option in newstart'
534              print*, ''
535#ifndef MESOINI
536              CALL ABORT
537#else
538              qsurf(ig,iq) = 0.
539#endif
540            endif
541             
542              if ((piceco2(ig) .eq. 0.).and.
543     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
544                     psolaralb(ig,1) = albedo_h2o_ice
545                     psolaralb(ig,2) = albedo_h2o_ice
546c                     PRINT*,'surfini.F frost',
547c     &                  lati(ig)*180./pi, long(ig)*180./pi
548               endif
549              enddo
550           endif
551          end do
552          PRINT*,'surfini: minimum albedo with frost and co2',
553     &            minval(albedodat)
554          PRINT*,'surfini: maximum albedo with frost and co2',
555     &            maxval(albedodat)
556       ELSE
557         watercaptag(:) = .false.
558       END IF ! OF IF(water)
559         
560      RETURN
561      END
Note: See TracBrowser for help on using the repository browser.