Changeset 3581 for trunk/LMDZ.MARS/libf/phymars/albedocaps.F90

Timestamp:

Jan 15, 2025, 6:35:44 PM (2 weeks ago)

Author:

jbclement

Message:

Mars PCM:

• Deletion of the variable 'albedo_co2_cap' which was redundant with 'albedo_perennialco2'.
• 'albedo_perennialco2' is now a vector of dimension 2 to account for different albedos of CO2 ice between the northern and southern hemispheres (by default 0.6 and 0.85 respectively).

JBC

File:

• trunk/LMDZ.MARS/libf/phymars/albedocaps.F90 (modified) (18 diffs)

trunk/LMDZ.MARS/libf/phymars/albedocaps.F90

@@ -11 +11 @@
 USE mod_phys_lmdz_transfert_para, ONLY: bcast
 USE mod_phys_lmdz_para, ONLY: is_master
-USE paleoclimate_mod, ONLY: paleoclimate,albedo_perennialco2,albedo_co2_cap
+USE paleoclimate_mod, ONLY: paleoclimate, albedo_perennialco2
 
 implicit none
@@ -42 +42 @@
 integer,parameter :: TESlonsize=72
 ! longitudes, in TES files, in degrees, from TESlon(1)=-177.5 to TESlon(72)=177.5
-real,save :: TESlon(TESlonsize) 
+real,save :: TESlon(TESlonsize)
 integer,parameter :: TESlatsize=30
 ! latitudes (north hemisphere file), in degrees, from TESlatn(1)=31,
@@ -67 +67 @@
 if (firstcall) then
   ! find out if user wants to use TES cap albedoes or not
-  TESicealbedo=.false. ! default value 
+  TESicealbedo=.false. ! default value
   write(*,*)" albedocaps: Use TES Cap albedoes ?"
   call getin_p("TESicealbedo",TESicealbedo)
   write(*,*)" albedocaps: TESicealbedo = ",TESicealbedo
-  
+
   ! if using TES albedoes, load coeffcients
   if (TESicealbedo) then
@@ -78 +78 @@
     call getin_p("TESice_Ncoef",TESice_Ncoef)
     write(*,*)" albedocaps: TESice_Ncoef = ",TESice_Ncoef
-    
+
     write(*,*)" albedocaps: Coefficient for Southern Cap ?"
     TESice_Scoef=1.0 ! default value
@@ -84 +84 @@
     write(*,*)" albedocaps: TESice_Scoef = ",TESice_Scoef
   endif
-  
+
   call read_TES_icecap_albedo(zls_old,pi,radeg,TESlon,TESlatn,TESlats,TESls,TESalbn,TESalbs)
-  
+
   firstcall=.false.
 endif ! of if (firstcall)
-    
-do ig=1,ngrid
-  if (latitude(ig).lt.0.) then
-    icap=2 ! Southern hemisphere
-  else
-    icap=1 ! Northern hemisphere
-  endif
-  
-  if (piceco2(ig).gt.0) then
-    ! set emissivity of surface to be the ice emissivity
-    emisref(ig)=emisice(icap)
-    ! set the surface albedo to be the ice albedo
-    if (TESicealbedo) then
-      call TES_icecap_albedo(zls,ig,psolaralb(ig,1),icap,zls_old,pi,radeg,TESlon,TESlatn,TESlats,TESls,TESalbn,TESalbs)
-      psolaralb(ig,2)=psolaralb(ig,1)
+
+do ig = 1,ngrid
+    if (latitude(ig) < 0.) then
+        icap = 2 ! Southern hemisphere
     else
-      if(paleoclimate) then
-            psolaralb(ig,1) = albedo_co2_cap
-            psolaralb(ig,2) = albedo_co2_cap
-      else
+        icap = 1 ! Northern hemisphere
+    endif
+
+    if (piceco2(ig) > 0) then ! CO2 frost
+        ! set emissivity of surface to be the ice emissivity
+        emisref(ig) = emisice(icap)
+        ! set the surface albedo to be the ice albedo
+        if (TESicealbedo) then
+            call TES_icecap_albedo(zls,ig,psolaralb(ig,1),icap,zls_old,pi,radeg,TESlon,TESlatn,TESlats,TESls,TESalbn,TESalbs)
+            psolaralb(ig,2) = psolaralb(ig,1)
+        else
             psolaralb(ig,1) = albedice(icap)
             psolaralb(ig,2) = albedice(icap)
-      endif
-    endif
-  else if(paleoclimate .and. piceco2_peren(ig) > 0.) then
-     psolaralb(ig,1) = albedo_perennialco2
-     psolaralb(ig,2) = albedo_perennialco2
-     emisref(ig)=emisice(icap)
-  else if (watercaptag(ig) .and. water) then
-    ! there is a water ice cap: set the surface albedo to the water ice one
-    ! to do : emissivity
-    emisref(ig) = 1
-    psolaralb(ig,1)=albedo_h2o_cap
-    psolaralb(ig,2)=albedo_h2o_cap
-  else
-    ! set emissivity of surface to be bare ground emissivity
-    emisref(ig)=emissiv
-    ! set the surface albedo to bare ground albedo 
-    psolaralb(ig,1)=albedodat(ig)
-    psolaralb(ig,2)=albedodat(ig)
-  endif ! of if (piceco2(ig).gt.0)
+        endif
+    else if (paleoclimate .and. piceco2_peren(ig) > 0.) then ! CO2 perennial ice
+        psolaralb(ig,1) = albedo_perennialco2(icap)
+        psolaralb(ig,2) = albedo_perennialco2(icap)
+        emisref(ig) = emisice(icap)
+    else if (watercaptag(ig) .and. water) then ! H2O ice in polar caps
+        ! there is a water ice cap: set the surface albedo to the water ice one
+        ! to do : emissivity
+        emisref(ig) = 1
+        psolaralb(ig,:) = albedo_h2o_cap
+    else ! Bare ground
+        ! set emissivity of surface to be bare ground emissivity
+        emisref(ig) = emissiv
+        ! set the surface albedo to bare ground albedo
+        psolaralb(ig,:) = albedodat(ig)
+    endif ! of if (piceco2(ig).gt.0)
 enddo ! of ig=1,ngrid
 
@@ -143 +136 @@
 USE mod_phys_lmdz_para, ONLY: is_master
 USE mod_phys_lmdz_transfert_para, ONLY: bcast
-                  
+
 implicit none
 
@@ -156 +149 @@
 integer,parameter:: TESlonsize=72
 ! longitudes, in TES files, in degrees, from TESlon(1)=-177.5 to TESlon(72)=177.5
-real:: TESlon(TESlonsize) 
+real:: TESlon(TESlonsize)
 integer,parameter:: TESlatsize=30
 ! latitudes (north hemisphere file), in degrees, from TESlatn(1)=31,
@@ -195 +188 @@
     write(*,*) "albedocaps: using file ",trim(datadir)//"/npsc_albedo.nc"
   ENDIF
-  
+
   ierr=nf90_inq_varid(nid,"longitude",nvarid)
   if (ierr.ne.NF90_NOERR) then
@@ -209 +202 @@
     endif
   endif
-  
+
   ierr=nf90_inq_varid(nid,"latitude",nvarid)
   if (ierr.ne.NF90_NOERR) then
@@ -223 +216 @@
     endif
   endif
-  
+
   ierr=nf90_inq_varid(nid,"time",nvarid)
   if (ierr.ne.NF90_NOERR) then
@@ -237 +230 @@
     endif
   endif
-  
+
   ierr=nf90_inq_varid(nid,"albedo",nvarid)
   if (ierr.ne.NF90_NOERR) then
@@ -304 +297 @@
 
   zls_old=-999 ! dummy initialization
-    
+
 endif !is_master
-    
+
     call bcast(TESlon)
     call bcast(TESlatn)
@@ -316 +309 @@
     call bcast(zls_old)
     call bcast(radeg)
-  
+
 end subroutine read_TES_icecap_albedo
 
@@ -327 +320 @@
                   nf90_strerror, nf90_inq_varid, nf90_get_var, nf90_close
 USE mod_phys_lmdz_transfert_para, ONLY: bcast
-                  
+
 implicit none
 
@@ -364 +357 @@
 integer,parameter:: TESlonsize=72
 ! longitudes, in TES files, in degrees, from TESlon(1)=-177.5 to TESlon(72)=177.5
-real:: TESlon(TESlonsize) 
+real:: TESlon(TESlonsize)
 integer,parameter:: TESlatsize=30
 ! latitudes (north hemisphere file), in degrees, from TESlatn(1)=31,
@@ -380 +373 @@
 real:: TESalbs(TESlonsize,TESlatsize,TESlssize)
 ! encompassing nodes arranged as follow : 4 3
-real :: val(4)                          ! 1 2 
+real :: val(4)                          ! 1 2
 
 ! 1. Identify encompassing latitudes
@@ -439 +432 @@
 if (zls.ne.zls_old) then
   zlsd=zls*radeg ! solar longitude, in degrees
-  
+
   if (zlsd.lt.TESls(1)) then
     tinf=TESlssize
@@ -461 +454 @@
     endif
   endif ! of if (zlsd.lt.TESls(1))
-  
+
   zls_old=zls ! store current zls
 endif ! of if (zls.ne.zls_old)